一階常微分方程式の初期値問題(前進オイラー法とリープフロッグ法の比較)

 \dfrac{dy}{dx}=-y-\pi \exp(-x)\sin\pi xを解く。
この方程式の解析解は、 e^{-x}\cos\pi x

clear all; close all;

dydx = @(x,y) -y -pi*exp(-x)*sin(pi*x);
N = 20;
x0 = 0;
xN = 2;

%%%% Forward Euler %%%%
x = linspace(x0,xN,N);
y = zeros(1,N);
h = (xN-x0)/N;
y(1) = 1;

for k=2:N   y(k) = y(k-1) + dydx(x(k-1),y(k-1))*h;
end
plot(x, y, 'DisplayName','Forward Euler');
hold on;
plot(x, exp(-x).*cos(pi*x), 'DisplayName','Exact Solution');

%%%% Leapfrog %%%%
x2 = linspace(x0,xN,2*N);
y2 = zeros(1,2*N);
h2 = (xN-x0)/(2*N);
y2(1) = 1;
y2(2) = y2(1) + h2*dydx(x2(1),y2(1));

for k=2:2*N-1
   y2(k+1) = y2(k-1) + 2*h2*dydx(x2(k),y2(k));
end
plot(x2(1:2:end), y2(1:2:end), 'DisplayName','LeapFrog');
[x3,y3]=ode45(@(x,y) -y-pi*exp(-x)*sin(pi*x),[0 2],1);
plot(x3,y3,'DisplayName','ode45');
legend();

刻み数N=20のとき
f:id:seinzumtode:20200518162246p:plain

刻み数N=40のとき
f:id:seinzumtode:20200518162241p:plain

Nを増加させていくと、Leapfrog法の方が厳密解に近づくことがわかる。