Matlab, slope field, euler ODE system and 2nd order Runge-Kutta

Question

I have system of differential equations

x' = ax - by

y' = bx + ay

I need to find approximate solution using explicit Euler method and second order Runge-Kutta, Conditions a = 0, b=1,

x(0) = 0, y(0) = 1, and accualy for more than that,

Using WolframAlpha, I determined exact solutions to this, and tried to implement algorithms on my own. However, im not sure what im doing wrong, because they seem to be too much dependent on h parameter and N. Below, I post my matlab code. I don't like the way my approximate solution is shifted or too sparse/dense compared to exact solution. I have to play with parameters N and h too much to get something decent, and im not sure if this is correct. This code should run right away if someone would like to help me. Thank you!

close all
clear all
[x,y] = meshgrid(-4:0.5:4);
a = 0;
% Task 9.a)
for b=[1 0.5 0.1 0.05]
dx = a*x-b*y;
dy = b*x+a*y;
figure();
quiver(x,y,dx,dy);
end

% task b,c for Explicit Euler's 
a = 0;
b = 1;

N = 1000; % 
h = 0.0125

t = linspace(0,4*pi,N);
xa = []
ya = []
xa(1) = 0;
ya(1) = 1;
    xe = cos(t) - sin(t); % exact solutions 
    ye = sin(t) + cos(t);

for h=[0.1 0.05 0.025 0.0125]    
    for k=1:N-1    
        xa(k+1) = xa(k) + h*(a*xa(k) - b*ya(k));
        ya(k+1) = ya(k) + h*(b*xa(k) + a*ya(k));   
    end

figure()
plot(xa);
hold on
plot(ya);
hold on
plot(xe,'r');
hold on
plot(ye,'b');

end


% runge-kutta
xa_rk = []
ya_rk = []
a = 0;
b = 1;

xa_rk(1) = 0
ya_rk(1) = 1
for h=[0.05 0.025 0.0125 0.005 0.0025 0.001]    
    for k=2:N-1    
        K1_1=  a*xa_rk(k-1) - b*ya_rk(k-1);
        K2_1 = b*xa_rk(k-1) + a*ya_rk(k-1);

        xintmid = xa_rk(k-1) + K1_1* h/2; % K2        
        yintmid = ya_rk(k-1) + K2_1* h/2; % K2

        K1_2 = a*xintmid - b*yintmid;
        K2_2 = b*xintmid + a*yintmid;

        xa_rk(k) = xa_rk(k-1) + h*(K1_1+K1_2)/2;
        ya_rk(k) = ya_rk(k-1) + h*(K2_1+K2_2)/2;

    end

figure()
plot(xa_rk);
hold on
plot(ya_rk);
hold on
plot(xe,'r');
hold on
plot(ye,'b');

end

Answer 1

You have two problems in your code:

The exact solution is incorrect

Using Wolfram Alpha, the exact solution is:

The plot domains are incorrect:

In the exact solution, your domain is [0,4*pi] with 1000 points in between.

In the numerical solutions, you are using different values of h, without changing the number of iterations. This means that the domain size is from 0 to h*N. This is incorrect. If you want to vary the value of h, you need to specify the number of iterations so that your last iteration is at the desired last point of the domain (4*pi, for example).

One way of doing this is to vary N, then compute h as:

h = 4*pi/N    # could be off-by-one ... I didn't check 

This produces matching results between the exact and numerical solutions.