Runge-Kutta 4th Order on a 2nd Order ODE
131 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Matthew Worker
am 8 Feb. 2021
Bearbeitet: Manoj Kumar Ghosh
am 22 Okt. 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 Kommentare
Akzeptierte Antwort
James Tursa
am 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 Kommentare
Weitere Antworten (1)
Manoj Kumar Ghosh
am 21 Okt. 2022
Bearbeitet: Manoj Kumar Ghosh
am 22 Okt. 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 Kommentare
Siehe auch
Kategorien
Mehr zu Ordinary Differential Equations finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!