% 常微分方程式の数値解法　実行サンプル
%   dy/dt = -cy （Dahlquistのテスト方程式）
clear all;
clc;
close all;

% MATLAB の出力（表示）桁数を倍精度にするおまじない
%format long;

% 諸パラメータ（ここをいろいろ変更してみるとよい）
ini = 1.0;      % 初期値（変えても大勢に影響ない）
c = -1.0;       % テスト方程式の係数（これとΔtの積が安定性を支配する）
et = 10.0;       % [0, et] を積分
n = 20;         % 分割数
dt = et / n;    % Δt = et/n
tv = linspace(0, et, n+1);      % t方向のグリッド

% 近似解の器と初期値を用意する
y = zeros(length(tv),1);        % 解の入れ物を用意しておく（MATLABでは実は必須ではない）
y(1) = ini;     % 初期値代入； y(1) に入るので添え字に注意すること

% 手法を選択
mtd = 1;           % 手法：1=陽的Euler，2=陰的Euler，3=中点則，4=陰的中点則

% メインループ
for i = 2:n+1       % 時刻 et の数値解は y(n+1) なので添え字の範囲に注意すること
    switch mtd
        case 1      % 陽的 Euler
            mtdtxt = '陽的 Euler';
            y(i) = y(i-1) + dt * c * y(i-1);
        case 2      % 陰的 Euler
            mtdtxt = '陰的 Euler';
            y(i) = y(i-1) / (1.0 - c * dt);
        case 3      % 中点則
            mtdtxt = '中点則';
            if (i ==2)      % 中点則は多段法なので「出発値」を別途作る必要がある
                y(i) = y(i-1) + dt * c * y(i-1);    % 陽的Eulerで作る
            else
                y(i) = y(i-2) + 2.0 * dt * c * y(i-1);
            end
        case 4      % 陰的中点則
            mtdtxt = '陰的中点則';       % 線形の問題では台形則と同じ
            y(i) = (1 + dt * c/2.0) / (1 - dt * c/2.0) * y(i-1);
    end
end

%% 以降，解の描画

% 厳密解を用意
tvf = linspace(0, et, 100);         % 厳密解のグリッド（細かめに刻む）
ey = ini * exp( c * tvf);           % 厳密解（解析解から計算）

% 安定領域に入るかどうかを見るべき数
z = c * dt;

% 描画
plot(tv, abs(y), 'r-o');     % 近似解
hold on;
plot(tvf, abs(ey), 'b-');    % 厳密解
legend(mtdtxt, '厳密解', 'FontSize', 14);
title([mtdtxt '; n = ' num2str(n) '; Δt = ' num2str(dt) '; cΔt = ' num2str(z)], 'FontSize', 14);
xlabel('t', 'FontSize', 14);
ylabel('y', 'FontSize', 14);
ylim([-0.5 1.5]);       % y 方向の描画範囲 


%% コメント
% 上のソースでは，問題や解法をメインソースに埋め込んでしまっており，
% 問題を変えようとすると全面的に書き換えることになる．
% 迷った末，このサンプルプログラムは（あまり複雑にしないため）そこまでしていないが，
% ある程度問題を色々試してみるのであれば，解法ルーチンをサブルーチンとして独立させ，
% 問題も関数として一般的に定義して，
% 任意の問題に任意の解法を試せるようなプログラムを書く方が良い．
% 自宅の部屋の整理・掃除と一緒で，一時期，ある程度の手間（時間）をかけても，
% その後生活しやすいように物事を最適化しておく方がゆくゆくは効率が良くなる．
% （ただし「人生でこの１回だけ」といったイベントは，ブルートフォースでその場で無理矢理
% 解いてしまって準備に手間をかけない方がお得であることが多い（※個人の感想です））．