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

Question

我正在为我的一个类做作业，我应该使用我选择的程序(我选择了Matlab)编写代码，以使用四阶龙格-库塔方法解决贝塞尔函数微分方程。参考贝塞尔函数DE为:

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

我可以把它分成两个耦合的一阶DEs:

(J_n) ' = Z_n和

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

在此作业之前，我没有Matlab或任何其他编程语言的经验。我知道Matlab有'ode45'命令，但我应该自己写代码，不依赖于Matlab的命令。到目前为止，我一直在研究贝塞尔函数n=0的情况，但是当我试图画出函数时，我总是得到一个错误。当前的错误我说:"未定义的函数或方法'J'的输入参数类型'double'。"但我不知道如何修复这个错误，也不知道我的代码是否正确。有人能告诉我哪里错了，或者写这段代码的正确方法是什么?

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。