bessel函数求导的问题
本帖最后由 allenfieldin 于 2010-10-29 09:38 编辑程序如下：
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()
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()
legend('L','dL')
title('besseli 函数及其导数')
grid on
作图如下：
图看不清楚
总感觉besseli函数图象画错了，是不是不能这样求他的导数？ 觉得是自变量不对，不是)/dr和/dr/xi 2# messenger
messenger，
besselk和besseli函数求解的步长分别是dr和xi*dr，不是吗？在差分后求导不对吗？ 对不对，找个例子求一下不就知道了？
函数y=x.^2>> f=inline('x.^2','x')
f =
Inline function:
f(x) = x.^21.步长为1>> f01=f()
f01 =
1 4 9 16 25 36 49 64
>> diff(f01)/1
ans =
3 5 7 9 11 13 152.步长0.1>> f02=f()
f02 =
1.0000 1.2100 1.4400 1.6900 1.9600 2.2500 2.5600 2.8900 3.2400
>> diff(f02)/.1
ans =
2.1000 2.3000 2.5000 2.7000 2.9000 3.1000 3.3000 3.50003.步长0.01>> f03=f()
f03 =
1.0000 1.0201 1.0404 1.0609 1.0816 1.1025 1.1236 1.1449 1.1664
>> diff(f03)/.01
ans =
2.0100 2.0300 2.0500 2.0700 2.0900 2.1100 2.1300 2.1500既然采用的是向前差分，那么得到的7个数据就是前7个源数据的导数，容易看出，明显不对，步长越大误差越大，导数定义中步长可是要趋近于零的哦。
建议用数值计算书中的三点微分或五点微分写一个程序再试试看。另外，好像薛定宇书中有一个这样的现成数值微分计算公式 4# bainhome
多谢！！ 4# bainhome
另一方面我也很怀疑，因为别人也用这样的方法求解，都可以啊，为什么这个besseli函数就不行呢？很怀疑，我还是看看这本书，再编一下吧 besseli的图形看起来没有错啊，你觉得那里不对呢。这样求导数是一种近似的方法，应该根据实际中具体问题具体分析。另外，bessel的函数有其特性，尤其是其导数，你不妨查查数学参考书，应该有提到的。 7# taohe
syms x
y=besseli(1,x);
y1=diff(y);
上面直接可以得到besseli函数的倒数表达式:
dy=besseli(0,x)-1/x*besseli(1,x)
而不需要那样每个步长差分求解
当然两种方法都做了，二者求得的结果相差是不大的。利用dy上面的这种方法当然更简单
页:
[1]