微分方程式の初期値問題（オイラー法での数値振動）

を考える。

clear; close all;

Ns=[8,10,15,20,30];
x0=0;
xN=20;

figure();
hold on;
for idx=1:length(Ns)
    N = Ns(idx);
    x=linspace(x0,xN,N);
    y=zeros(1,N);
    y(1)=0.1;
    h=(xN-x0)/N;
    f = @(y) y-y^2;
    for k=2:N
        y(k) = y(k-1) + h*(f(y(k-1)));
    end    
    plot(1:N,y,'DisplayName',sprintf('N=%d',N));
    xlim([0,10]);
end
[x2,y2]=ode45(@(x,y)y-y^2,[0 20],0.1);
plot(x2,y2,'DisplayName','ode45');
legend();

N=8：不規則な振動
N=10：規則的な振動
N=15,20,30：振動なし

ロジスティック写像を用いた振動の説明
とおくと、が成り立つ

clear; close all;

x0=0;
xN=20;
Ns=[6,8,10,15,20,30];

figure(1);
for idx=1:length(Ns)
    N = Ns(idx);
    h=(xN-x0)/N;
    alpha = 1+h;
    
    x=linspace(x0,xN,N);
    y=zeros(1,N);
    z=zeros(1,N);
    y(1)=0.1;
    h=(xN-x0)/N;
    f = @(y) y-y^2;
    for k=2:N
        y(k) = y(k-1) + h*(f(y(k-1)));
    end    
    for k=1:N
        z(k)=h*y(k)/(1+h);
    end
    
    subplot(2,3,idx);
    hold on;
    plot(z(1:end-1),z(2:end),'ro');
    for k=1:N-1
        txt = strcat('\leftarrow',sprintf(' z(%d)',k));
        text(z(k),z(k+1),txt);
    end
    
    m = polyfit(z(1:end-1),z(2:end),2);
    func = @(x) m(1)*x.^2+m(2)*x+m(3);
    t=0:0.1:1;
    plot(t,func(t),'b');
    plot(0:1,0:1,'k-');
    title(sprintf('N=%d',N));    
end

figure(1);
sgtitle('ロジスティック写像を用いた区分幅hの大小による微分方程式の振動の説明');

N≥15では、直線y=xと曲線の交点（不動点）の周りに解が収束している。