Octaveでのオイラー法とルンゲクッタ法について
講義ノートを見ながら数値解析ソフトOctaveの操作を勉強しています。
読んでいる講義ノートにオイラー法とルンゲ・クッタ法でf''(x)+f'(x)+x=0の微分方程式を解く箇所があるのですが、書いている通りに打ち込むとx=0の直線になってしまいます。lsodeで解くとうまくいくのですが。学び始めたばかりの初心者でまったく分からないので、、どうすればいいか教えてください。以下が講義ノートに書いてあったプログラムです。xdot2,euler,rk4は定義された関数で、rk4eulerがメインプログラムです。
% xdot2.m
function xdot2 = xdot2(x,t)
xdot2=[x(2) -x(2)-x(1)];
% euler.m
function x=euler(F,x,t)
[m n]=size(t); %t のサイズ(=ステップ数)を抽出
for k=1:m-1 %tのサイズ-1回繰り返す。
tau=t(k+1)-t(k);
x(k+1,:)=x(k,:)+tau.*feval(F,x(k,:),t(k));
%上の式がオイラー法の公式
end
% rk4.m
function x=rk4(F,x,t)
[m n]=size(t);
for k=1:m-1 %t のサイズ-1 だけステップを繰り返す。
tau = t(k+1)-t(k); %刻み幅
F1=feval(F,x(k,:),t(k));
F2=feval(F,x(k,:)+1/2*tau.*F1,t(k)+1/2*tau);
F3=feval(F,x(k,:)+1/2*tau.*F2,t(k)+1/2*tau);
F4=feval(F,x(k,:)+tau*F3,t(k)+tau);
x(k+1,:)= x(k,:)+tau/6*(F1 +2*F2 +2*F3 +F4);
%以上が4次のRK法の公式
end
%rk4euler.m %オイラー、RK、lsode の比較
t=[0:0.2:10]’; x0=[0 1]; %定義域、初期値の定義
R1=euler('xdot2',x0,t); x1=R1(:,1); %オイラー法
R2= rk4('xdot2',x0,t); x2=R2(:,1); %RK 法
R3=lsode('xdot2',x0,t); x3=R3(:,1); %lsode
plot(t,x1,'rx-.',t,x2,'bx',t,x3,'ko');
xlabel('t'); ylabel('x');
legend('Euler','RK4','lsode');
title('rk4euler1.m');
