bessel函数求导的问题

allenfieldin

程序如下：
clear all
clc
r=100;
dr=0.1;
xi=0.7;
L=120;
Nr=r/dr;
r_1=(1:2*Nr)*dr;          %besselk函数的求解域
r_2=r_1*xi;                   %besseli函数的求解域
K=besselk(1,r_1);       %bessel函数的阶数为1
L=besseli(1,r_2);          %besseli函数的阶数为1
dK=diff(K)/dr;                %besselk函数的求导
dL=diff(L)/dr/xi;            %besseli函数的求导
subplot(1,2,1)
plot(r_1,K,'-r')
hold on
plot(r_1(1:2*Nr-1),dK,'.-r',r_1,0,'-k')
axis([0 4 -5 5])
legend('K','dB')
title('besselk函数及其导数')
grid on
subplot(1,2,2)
plot(r_2,L,'-b')
hold on
plot(r_2(1:2*Nr-1),dL,'.-b',r_2,0,'-k')
axis([0 4  -2 10])
legend('L','dL')
title('besseli 函数及其导数')
grid on

作图如下：
图看不清楚
总感觉besseli函数图象画错了，是不是不能这样求他的导数？

messenger

觉得是自变量不对，不是)/dr和/dr/xi

allenfieldin

2# messenger
messenger，
besselk和besseli函数求解的步长分别是dr和xi*dr，不是吗？在差分后求导不对吗？

bainhome

对不对，找个例子求一下不就知道了？
函数y=x.^2

1. >> f=inline('x.^2','x')
   
2. f =
   
3.      Inline function:
   
4.      f(x) = x.^2

复制代码

1.步长为1

1. >> f01=f([1:8])
   
2. f01 =
   
3.      1     4     9    16    25    36    49    64
   
4. >> diff(f01)/1
   
5. ans =
   
6.      3     5     7     9    11    13    15

复制代码

2.步长0.1

1. >> f02=f([1:.1:1.8])
   
2. f02 =
   
3.     1.0000    1.2100    1.4400    1.6900    1.9600    2.2500    2.5600    2.8900    3.2400
   
4. >> diff(f02)/.1
   
5. ans =
   
6.     2.1000    2.3000    2.5000    2.7000    2.9000    3.1000    3.3000    3.5000

复制代码

3.步长0.01

1. >> f03=f([1:.01:1.08])
   
2. f03 =
   
3.     1.0000    1.0201    1.0404    1.0609    1.0816    1.1025    1.1236    1.1449    1.1664
   
4. >> diff(f03)/.01
   
5. ans =
   
6.     2.0100    2.0300    2.0500    2.0700    2.0900    2.1100    2.1300    2.1500

复制代码

既然采用的是向前差分，那么得到的7个数据就是前7个源数据的导数，容易看出，明显不对，步长越大误差越大，导数定义中步长可是要趋近于零的哦。
建议用数值计算书中的三点微分或五点微分写一个程序再试试看。另外，好像薛定宇书中有一个这样的现成数值微分计算公式

allenfieldin

4# bainhome
多谢！！

allenfieldin

4# bainhome
另一方面我也很怀疑，因为别人也用这样的方法求解，都可以啊，为什么这个besseli函数就不行呢？很怀疑，我还是看看这本书，再编一下吧

taohe

besseli的图形看起来没有错啊，你觉得那里不对呢。这样求导数是一种近似的方法，应该根据实际中具体问题具体分析。另外，bessel的函数有其特性，尤其是其导数，你不妨查查数学参考书，应该有提到的。

allenfieldin

7# taohe
syms x
y=besseli(1,x);
y1=diff(y);

上面直接可以得到besseli函数的倒数表达式:
dy=besseli(0,x)-1/x*besseli(1,x)
而不需要那样每个步长差分求解
当然两种方法都做了，二者求得的结果相差是不大的。利用dy上面的这种方法当然更简单