How to plot a phase portrait using Runge-Kutta for a system of ODEs?

Question

Iris 2024-3-24

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2098071-how-to-plot-a-phase-portrait-using-runge-kutta-for-a-system-of-odes

评论： Iris 2024-3-24

采纳的回答： Sam Chak

在 MATLAB Online 中打开

Hello, I would like to plot a phase portrait for the following system of ODEs:

y'=z

z'=y(1-y^2)

I want to do this manually, without using ode45 or some other built in solver, so I'm trying to build out a Runge-Kutta solver that I can use to make the phase portrait. I did this problem by hand, and the plot I am getting from my code does not match the expected answer. Can anyone point me towards where I have gone wrong? Here is my code for the RK solver:

% input arguments: y',z', c (vector of initial values), t (linspace vector
% equally spaced points on [a,b])
n = length(t)-1; m = length(c); dt = t(2)-t(1);
w = zeros(m,n+1); w(:,1) = c; 
for i = 1:n
    k1 = dt*(w(:,i));
    k2 = dt*(w(:,i)+(k1/2));
    k3 = dt*(w(:,i)+(k2/2));
    k4 = dt*(w(:,i)+k3);
    k = (k1+(2*k2)+(2*k3)+k4)/6;
    w(:,i+1) = w(:,i) + k;
end

And for the phase portrait:

f = @(y,z) [z,y*(1-y^2)]; 
t = linspace(0,10,50);
hold on 
for y0 = -5:1:5 
    for z0 = -5:1:5
        c = [y0,z0]; % initial conditions 
        w = RK_solver(f,c,t); 
        y = w(1,:);
        z = w(2,:);
        plot(y,z)
        axis([-5 5 -5 5]);
    end
end

This is the picture I get, but there shoud be two centers and one saddle if I'm not mistaken- I don't know why this is all linear.

I am a bit rusty on differential equations. Any help is appreciated.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Sam Chak 2024-3-24

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2098071-how-to-plot-a-phase-portrait-using-runge-kutta-for-a-system-of-odes#answer_1430066

编辑：Sam Chak 2024-3-24

在 MATLAB Online 中打开

Hi @Iris

Are you expecting the followng Phase Portraits?

odefcn = @(t, y) [y(2); % y' = z

y(1)*(1 - y(1)^2)]; % z' = y·(1 - y²)

tspan = [0, 10];

y0 = -0.6:0.15:0.6;

z0 = y0;

for i = 1:numel(y0)

for j = 1:numel(z0)

[t, y] = ode45(odefcn, tspan, [y0(i); z0(j)]);

plot(y(:,1), y(:,2)), hold on

end

grid on, axis([-1.6 1.6 -1.6 1.6])

xlabel('y'), ylabel('z'), title('Phase Portraits')

5 个评论
显示 3更早的评论隐藏 3更早的评论

Sam Chak 2024-3-24

You're welcome, @Iris 👍. I'm just curious, is this particular topic part of a math course that focuses on solving differential Equations, or is it a more general MATLAB course that deals with various computational problems?

Iris 2024-3-24