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

2 次查看(过去 30 天)
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 个评论
James Tursa
James Tursa 2022-11-17
What have you done so far? What specific problems are you having with your code?
Cubii4
Cubii4 2022-11-27
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.

请先登录,再进行评论。

采纳的回答

Torsten
Torsten 2022-11-27
移动:Cris LaPierre 2022-11-28
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');

更多回答(1 个)

Cris LaPierre
Cris LaPierre 2022-11-10
  1 个评论
Cubii4
Cubii4 2022-11-27
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.

请先登录,再进行评论。

产品


版本

R2015a

Community Treasure Hunt

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

Start Hunting!

Translated by