Need to solve this with Runge-Kutta 4th order method in MatLab.

2 vues (au cours des 30 derniers jours)
Cubii4
Cubii4 le 10 Nov 2022
y''-μ(1-y^2)*y'+y=0
y(x=0)=1 ; y'(x=0)=0
μ=0; 0.2; 1; 5 and 10 (for each value a different result)
y'(x) versus y(x) and x versus y(x) functions plot in the range [0; 50].
  2 commentaires
James Tursa
James Tursa le 17 Nov 2022
What have you done so far? What specific problems are you having with your code?
Cubii4
Cubii4 le 27 Nov 2022
fy=@(x,y,z) z;
fz=@(x,y,z) (1-y^2)*z-y;
x(1)=0;
z(1)=0;
y(1)=1;
h=1;
xfinal=10;
N=ceil((xfinal-x(1))/h);
for j=1:N:xfinal
x(j+h)=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+h)=y(j)+h/6*(k1y+2*k2y+2*k3y+k4y);
z(j+h)=z(j)+h/6*(k1z+2*k2z+2*k3z+k4z);
end
x = 1:N:xfinal;
figure;
plot(x,y,z,'LineWidth',2);
xlabel('X');
ylabel('Y');
I took out the μ value and decreased the range to (0,10) just until I can understand how to form the code and this is what I have done but there is a problem with the plot. Can I get some help with it. Thanks.

Connectez-vous pour commenter.

Réponse acceptée

Torsten
Torsten le 27 Nov 2022
Déplacé(e) : Cris LaPierre le 28 Nov 2022
fy=@(x,y,z) z;
fz=@(x,y,z) (1-y^2)*z-y;
x(1)=0;
z(1)=0;
y(1)=1;
h=0.1;
xfinal=10;
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
figure;
plot(x,[y;z],'LineWidth',2);
xlabel('X');
ylabel('Y');

Plus de réponses (1)

Cris LaPierre
Cris LaPierre le 10 Nov 2022
  1 commentaire
Cubii4
Cubii4 le 27 Nov 2022
I still have some problems with the plot. I wrote the code in the comment above if I can get some help with it. Thanks.

Connectez-vous pour commenter.

Produits


Version

R2015a

Community Treasure Hunt

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

Start Hunting!

Translated by