Решение функции Бесселя с помощью Runge Kutta

Question

Я работаю над заданием для моего класса, и я должен написать код, используя программу по своему выбору (я выбрал Matlab) для решения дифференциального уравнения функции Бесселя используя метод Рунге-Кутты 4-го порядка. Для справки функция Бесселя DE:

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

Я способен отделить это в двух соединенных первого порядка КЭУ по:

(J_n) '= Z_n и

(Z_n) + (1/х) * Z_n + [(х^2-п^2)/х^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?)

Кроме того, не забудьте лет ur index i при вычислении dJ и dZ в цикле for.