Matlab、斜率场、欧拉 ODE 系统和二阶 Runge-Kutta

Question

我有微分方程组

x' = ax - by

y' = bx + ay

我需要使用显式欧拉法和二阶龙格-库塔法找到近似解，条件 a = 0, b=1,

x(0) = 0, y(0) = 1，而且实际上不止于此，

使用 WolframAlpha，我确定了对此的精确解决方案，并尝试自己实现算法。但是，我不确定我做错了什么，因为它们似乎过于依赖 h 参数和 N。下面，我 post 我的 matlab 代码。与精确解相比，我不喜欢我的近似解偏移的方式或 sparse/dense。我必须过多地使用参数 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，精确解是：

$x=-sin(t) ,\qquad y=cos(t)$

绘图域不正确：

在精确解中，您的域是 [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

这会在精确解和数值解之间产生匹配结果。

Matlab、斜率场、欧拉 ODE 系统和二阶 Runge-Kutta

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

matlab

system

ode

differential-equations

runge-kutta

精确解不正确

绘图域不正确：