2階常微分方程式の数値解法

を解く。
初期条件はおよび
上記は以下の連立微分方程式に変形できる。

clear all; close all;

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

%%%% Forward Euler %%%%
x = linspace(x0,xN,2*N);
y = zeros(1,2*N);
z = zeros(1,2*N);
h = (xN-x0)/(2*N);
y(1) = 1;
z(1) = 0;
for k=2:2*N-1
    y(k) = y(k-1) + h*z(k-1);
    z(k) = z(k-1) - h*y(k-1);
end
plot(x, y, 'DisplayName','Forward Euler');
hold on;
plot(x, cos(x), 'DisplayName','Exact Solution');

%%%% Leapfrog %%%%
x2 = linspace(x0,xN,2*N);
y2 = zeros(1,2*N);
z2 = zeros(1,2*N);
h2 = (xN-x0)/(2*N);
y2(1) = 1;
z2(1) = 0;
y2(2) = y2(1)+h*z2(1);
z2(2) = z2(1)-h*y2(1);
for k=1:N-1
   y2(2*k+1) = y2(2*k-1)+2*h*z2(2*k);
   if(k<N-1)
     z2(2*k+2) = z2(2*k)-2*h*y2(2*k+1);
   end
end
plot(x2(1:2:end), y2(1:2:end), 'DisplayName','LeapFrog');

%%%% matlab solver %%%%
syms y(x);
Dy = diff(y);
ode = diff(y,x,2) == -y;
cond1 = y(0) == 1;
cond2 = Dy(0) == 0;
conds = [cond1 cond2];
x3=linspace(x0,xN,2*N);
y3(x) = dsolve(ode,conds);
plot(x3,y3(x3),'DisplayName','dsolve');
legend();

title('2nd Order ODE');

位相図を描いてみる。

clear all; close all;

N = 40;
x0 = 0;
xN = 10;

%%%% Exact Solution %%%%
x = linspace(x0,xN,2*N);
y = cos(x);
z = -sin(x);
plot(y,z,'DisplayName','Exact Solution');
hold on;
axis equal;

%%%% Forward Euler %%%%
x = linspace(x0,xN,2*N);
y = zeros(1,2*N);
z = zeros(1,2*N);
h = (xN-x0)/(2*N);
y(1) = 1;
z(1) = 0;
for k=2:2*N-1
    y(k) = y(k-1) + h*z(k-1);
    z(k) = z(k-1) - h*y(k-1);
end
plot(y(1:end-1),z(1:end-1),'DisplayName','Forward Euler');

%%%% Leapfrog %%%%
x2 = linspace(x0,xN,2*N);
y2 = zeros(1,2*N);
z2 = zeros(1,2*N);
h2 = (xN-x0)/(2*N);
y2(1) = 1;
z2(1) = 0;
y2(2) = y2(1)+h*z2(1);
z2(2) = z2(1)-h*y2(1);
for k=1:N-1
   y2(2*k+1) = y2(2*k-1)+2*h*z2(2*k);
   if(k<N-1)
     z2(2*k+2) = z2(2*k)-2*h*y2(2*k+1);
   end
end
plot(y2(1:2:end-1),z2(2:2:end),'DisplayName','LeapFrog');

legend();