Runge-Kutta 4th Order on a 2nd Order ODE

99 views (last 30 days)
I am dealing with the following 2nd order ODE:
With initial conditions of: , , .
and should be 3.24
I have found the system as:
For the life of me I cannot see to get the below code to produce the correct output... is it the code or my calculus?
clear all;
fy=@(x,y,z) 6*x*z-5*z
fz=@(x,y,z) z;
x(1)=0;
z(1)=2/3;
y(1)=0;
h=0.5;
xfinal=3;
N=ceil((xfinal-x(1))/h);
for j=1:N
x(j+1)=x(j)+h;
k1y=fy(x(j),y(j),z(j));
k1z=fz(x(j),y(j),z(j));
k2y=fy(x(j)+h/2,y(j)+h/2*k1y,z(j)+h/2*k1z);
k2z=fz(x(j)+h/2,y(j)+h/2*k1y,z(j)+h/2*k1z);
k3y=fy(x(j)+h/2,y(j)+h/2*k2y,z(j)+h/2*k2z);
k3z=fz(x(j)+h/2,y(j)+h/2*k2y,z(j)+h/2*k2z);
k4y=fy(x(j)+h,y(j)+h*k3y,z(j)+h*k3z);
k4z=fz(x(j)+h,y(j)+h*k3y,z(j)+h*k3z);
y(j+1)=y(j)+h/6*(k1y+2*k2y+2*k3y+k4y);
z(j+1)=z(j)+h/6*(k1z+2*k2z+2*k3z+k4z);
end
disp(y(N));
figure;
plot(x,y,'LineWidth',2);
xlabel('X');
ylabel('Y');

Accepted Answer

James Tursa
James Tursa on 8 Feb 2021
All of your derivatives should be with respect to the independent variable x. In the context of your problem, dy/dz does not make any sense. So you have two variables y and z, and two derivatives dy/dx and dz/dx. So your two 1st order equations are:
dy/dx = y' = z
dz/dx = d(y')/dx = y'' = (5*x*y' - 5*y) / (3*x^2) = (5*x*z - 5*y) / (3*x^2)
The dz/dx expression is obtained by setting your differential equation to 0 (you didn't specify this but I assumed it to be the case) and then solving for y'' and substituting in z for y'. So your function handles should be:
fy=@(x,y,z) z
fz=@(x,y,z) (5*x*z - 5*y) / (3*x^2)

More Answers (1)

Manoj Kumar Ghosh
Manoj Kumar Ghosh on 21 Oct 2022
Edited: Manoj Kumar Ghosh on 22 Oct 2022
clc
clear all
n=10000;
x=linspace(1,3,n+1);
y0=[0;2/3];
h=(x(end)-x(1))/n;
f=@(x,y,z) ([z;((5*z)/(3*x)-(5*y)/(3*x.^2))]) %took y'=z
for i=1:n
k1=h*f(x(i),y0(1),y0(2));
k2=h*f(x(i)+h/2,(y0(1)+k1(1)/2),(y0(2)+k1(2)/2));
k3=h*f(x(i)+h/2,(y0(1)+k2(1)/2),(y0(2)+k2(2)/2));
k4=h*f(x(i)+h,(y0(1)+k3(1)),(y0(2)+k3(2)));
y0=y0+(1/6)*(k1+2*k2+2*k3+k4);
end
disp(y0(1))

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!