Runge-Kutta 4th Order on a 2nd Order ODE
131 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Matthew Worker
le 8 Fév 2021
Modifié(e) : Manoj Kumar Ghosh
le 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 commentaires
Réponse acceptée
James Tursa
le 8 Fév 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 commentaires
Plus de réponses (1)
Manoj Kumar Ghosh
le 21 Oct 2022
Modifié(e) : Manoj Kumar Ghosh
le 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 commentaires
Voir également
Catégories
En savoir plus sur Ordinary Differential Equations dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!