4次のルンゲクッタ法による常微分方程式の数値解法

を計算する。
解析解は,

rungekutta.m

clear; close all;

x0=0;
xN=1;
N=10;
h=(xN-x0)/N;
x=linspace(x0,xN,N);
y=zeros(1,N);
a = 1;
y(1)=1;
f = @(x,y) a*y;

for k=2:N
  k1 = h*f(x(k-1),y(k-1));
  k2 = h*f(x(k-1)+h/2,y(k-1)+k1/2);
  k3 = h*f(x(k-1)+h/2,y(k-1)+k2/2);
  k4 = h*f(x(k-1)+h,y(k-1)+k3);
  y(k) = y(k-1) + 1/6*(k1+2*k2+2*k3+k4);
end
plot(x,y,'DisplayName','Runge Kutta');
hold on;
exact = exp(a*x);
plot(x,exact,'DisplayName','Exact Solution');

y2=zeros(1,N);
y2(1)=1;
for k=2:N
  k1 = h*f(x(k-1),y2(k-1));
  k2 = h*f(x(k-1)+h,y2(k-1)+k1);
  y2(k) = y2(k-1) + 1/2*(k1+k2);
end
plot(x,y2,'DisplayName','Heun method');

legend();

fprintf('Runge Kutta RMS=%d\n',rms(y-exact));
fprintf('Heun Method RMS=%d\n',rms(y2-exact));

ホイン法（2次のルンゲクッタ）とあんまり変わらない

>> rungekutta
Runge Kutta RMS=1.298050e-01
Heun Method RMS=1.315339e-01

次は、
を計算する。
解析解は、

clear; close all;

x0=0;
xN=10;
N=10;
h=(xN-x0)/N;
x=linspace(x0,xN,N);
y=zeros(1,N);
a = 1/4;
b = pi/6;
c = 0.6;
y(1)=0.3;
f = @(x,y) a*y+c*exp(a*x)*cos(x+b);

for k=2:N
  k1 = h*f(x(k-1),y(k-1));
  k2 = h*f(x(k-1)+h/2,y(k-1)+k1/2);
  k3 = h*f(x(k-1)+h/2,y(k-1)+k2/2);
  k4 = h*f(x(k-1)+h,y(k-1)+k3);
  y(k) = y(k-1) + 1/6*(k1+2*k2+2*k3+k4);
end
plot(x,y,'DisplayName','Runge Kutta');
hold on;

y2=zeros(1,N);
y2(1)=1;
for k=2:N
  k1 = h*f(x(k-1),y2(k-1));
  k2 = h*f(x(k-1)+h,y2(k-1)+k1);
  y2(k) = y2(k-1) + 1/2*(k1+k2);
end
plot(x,y2,'DisplayName','Heun method');

%%% Matlab solver
syms y(t);
ode = diff(y,t) == a*y+c*exp(a*t)*cos(t+b);
cond = y(0) == 0.3;
ySol(t) = dsolve(ode,cond);
plot(x,eval(ySol(x)),'DisplayName','Matlab solver');

%%% Exact solution
plot(x,0.6*exp(x/4).*sin(x+pi/6),'DisplayName','Exact Solution');

legend();

ルンゲクッタが精度良く計算できていることがわかる。