使用龙格库塔求解贝塞尔函数

Question

我正在为我的一个班级做作业，我应该使用我选择的程序（我选择 Matlab）编写代码，以使用 4 阶龙格求解贝塞尔函数微分方程-库塔法。作为参考，贝塞尔函数 DE 是：

x^2*(J_n)''+x*(J_n)'+(x^2-n^2)*J_n=0。

我可以通过以下方式将其分成两个耦合的一阶 DE：

(J_n)'=Z_n 和

(Z_n)'+(1/x)*Z_n+[(x^2-n^2)/x^2]*J_n=0。

在此任务之前，我没有使用 Matlab 或任何其他编程语言的经验。我知道 Matlab 有“ode45”命令，但我应该自己编写代码，而不是依赖于 Matlab 的命令。到目前为止，我一直在研究 Bessel 函数的 n=0 情况，但是当我尝试绘制函数时，我一直收到错误消息。我说的当前错误是：“'double' 类型的输入参数的未定义函数或方法 'J'。”但我不知道如何解决这个错误，也不知道我的代码是否正确。有人可以告诉我哪里出了问题或者编写此代码的正确方法是什么？

h=0.01;             %step size

J_0(1)=1;           %initial condition for J_0

Z_0(1)=1;           %initial condition for Z_0-This value should be zero
                    %but Matlab gives me an error. To fix this, I input
                    %Z_0(1)-1 to use the correct value for Z_0(1).

x(1)=0.001;         %first value of x

dZ(Z_0,J_0)=(-1/x)*(Z_0-1)-J_0;

for i=[1:1:10]

    dZ1=(-1/x)*(Z_0-1)-J_0;
    dJ1=(Z_0(1)-1)*h;

    dZ2=(-1/x)*(Z_0-1+0.5*h)-(J_0+0.5*h*dJ1);
    dJ2=((Z_0(1)-1)+dZ1)*h;

    dZ3=(-1/x)*(Z_0-1+0.5*h)-(J_0+0.5*h*dJ2);
    dJ3=((Z_0(1)-1)+dZ1+dZ2)*h;

    dZ4=(-1/x)*(Z_0-1+h)-(J_0+h*dJ3);
    dJ4=((Z_0(1)-1)+dZ1+dZ2+dZ3)*h;

    J(i+1)=J(i)+(h/6)*(dJ1+2*dJ2+2*dJ3+dJ4);
end

plot(J_0);

提前感谢您的帮助

Answer 1

你的问题就出来了：

J(i+1)=J(i)+(h/6)*(dJ1+2*dJ2+2*dJ3+dJ4);

在赋值运算符的右侧，您使用变量 J 在 i 取值 1 之前从未设置过。对我来说看起来像是一个错字（应该改为 J_0 吗？）

另外，在 for 循环中计算 dJ 和 dZ 内容时，不要忘记索引 i。