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

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.

Answer 1

[Sam Chak]

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
end
grid on, axis([-1.6 1.6 -1.6 1.6])
xlabel('y'), ylabel('z'), title('Phase Portraits')

[Iris]

I think so. This is the (very) rough sketch I did, they look fairly similar

[Sam Chak]

Hi @Iris, The initial values of y0 and z0 spread out too wide. I updated the code to demo this.

[Iris]

Oh I see, perfect. I believe I fixed my Runge Kutta code and was able to achieve a similar result. Thank you very much for your demo.

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.*f(t(i),w(1,i),w(2,i));
    k2 = dt.*f(t(i)+dt/2,w(1,i)+(k1(1)/2),w(2,i)+(k1(2)/2));
    k3 = dt.*f(t(i)+dt/2,w(1,i)+(k2(1)/2),w(2,i)+(k2(2)/2));
    k4 = dt.*f(t(i)+dt,w(1,i)+k3(1),w(2,i)+k3(2));
    k = (k1+(2*k2)+(2*k3)+k4)/6;
    w(:,i+1) = w(:,i) + k;
end

.

[Sam Chak]

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]

The second- this is for a numerical analysis course.