Runge-Kutta 4th Order on a 2nd Order ODE
24 views (last 30 days)
Show older comments
Matthew Worker
on 8 Feb 2021
Edited: Manoj Kumar Ghosh
on 22 Oct 2022
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');
0 Comments
Accepted Answer
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)
0 Comments
More Answers (1)
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))
0 Comments
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!