Matlab, поле склона, система ODE Эйлера и Рунге-Кутта 2-го порядка

Question

У меня есть система дифференциальных уравнений

х ' = топор - на

y' = bx + ay

Мне нужно найти приближенное решение, используя явный метод Эйлера и Рунге-Кутта второго порядка, условия a = 0, b = 1,

x (0) = 0, y (0) = 1, и, более того,

Используя WolframAlpha, я определил точные решения для этого и попытался реализовать алгоритмы самостоятельно. Тем не менее, я не уверен, что я делаю неправильно, потому что они, кажется, слишком сильно зависят от параметра h и N. Ниже я публикую свой код matlab. Мне не нравится, как мое приближенное решение смещено или слишком разреженно / плотно по сравнению с точным решением. Я должен играть с параметрами N и h слишком много, чтобы получить что-то приличное, и я не уверен, что это правильно. Этот код должен запускаться сразу, если кто-то захочет мне помочь. Спасибо!

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

У вас есть две проблемы в вашем коде:

Точное решение неверно

Используя Wolfram Alpha, точное решение:

Домены сюжета неверны:

В точном решении, ваш домен [0,4*pi] с 1000 баллов между ними.

В численных решениях вы используете разные значения h без изменения количества итераций. Это означает, что размер домена от 0 до h*N, Это неверно Если вы хотите изменить значение h вам нужно указать количество итераций, чтобы ваша последняя итерация была в нужной последней точке домена (4*pi, например).

Один из способов сделать это - изменить N, а затем вычислить h как:

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

Это приводит к совпадению результатов между точным и численным решениями.