• 締切済み

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');

みんなの回答

  • f272
  • ベストアンサー率46% (8024/17152)
回答No.1

おかしなところは rk4euler.mの2行目で t=[0:0.2:10]’; x0=[0 1]; %定義域、初期値の定義 に変な文字「’」がある。 これを「'」に変えれば,ちゃんと動く。

関連するQ&A

  • オイラー法

    区間[0,10]においてh=0.01としオイラー法を用いてとけ。また誤差も出せ。 y'=3x^2-x^3-x^6+(2x^3+1)y-y^2 y(0)=0.5 なんですけどa=3くらいから値がおかしくなりa=10の時には誤差ですぎます。誤差がでないようにするにはどこを直せばいいのか教えてください。 おねがいします。 #include<stdio.h> #include<math.h> double h=0.01; double x0=0.5; main(void) { int a; double euler(int a); double gosa; for(a=1;a<=10;a++) gosa=shin(a)-euler(a); printf("%d %f %f\n",a,euler(a),gosa); } double shin(double a) { double y; y=a*a*a+1.0/(1.0+exp(-a)); return y; } double f(double t,double x){ double y; y = 3.0*t*t-t*t*t-t*t*t*t*t*t+(2.0*t*t*t+1.0)*x-x*x; return y; } /* Euler Function */ double euler(int a){// x0:Initial Value double T,nx,x; // Initial Value x = x0; // Euler Method for(T = 0.0;T < a;T += h){ nx = x + h*f(T,x); x = nx; }

  • オイラー法、ルンゲクッタ法について。

    オイラー法、ルンゲクッタ法について。 この2つについて分からない事があるので質問します。 まず、オイラーについてですが、yi+1=yi+hf(x,y)という式がテイラー展開によって求まると言われましたが、テイラー展開の2次以降の項は微少量として無視できるのは分かります。でもそもそもテイラー展開ってひとつ先の値を今の値から求まるみたいな展開でしたっけ??というのが一つ目の質問です。 2つ目は、オイラーの式の中のf(x,y)についてです。簡単なバネ・マス・ダンパ系を考えた時、運動方程式はm・d2x/dt2+c・dx/dt+kx=0となると思いますが、この場合のf(x,y)はどうやって求めるのでしょうか。 3つ目はルンゲクッタそもそもについてです。 ルンゲクッタとはK1K2K3K4という係数(?)に1221という重みをかけるとyi+1が求まるそうですが、この理由がどんなサイトや本を見ても納得出来ません。 何か分かりやすい本やサイトがあれば教えて頂けないでしょうか。 以上3つの質問、回答よろしくお願いします。

  • オイラー法、2次ルンゲクッタ法、4次ルンゲクッタ法のC言語プログラムに

    オイラー法、2次ルンゲクッタ法、4次ルンゲクッタ法のC言語プログラムについて教えてください! 課題なのですが、まったくわからず困ってます>< 1 常微分方程式 dy/dx=f(x,y),y(0)=1 の数値解をオイラー法を用いて計算するプログラムを作為せよ。ただし、f(x,y)=3-6x^2-4x+2xyとする。 2 α=1,β=1,γ=1/2,σ=1/2 の場合の2次ルンゲクッタ法を考える。1と同じ常微分方程式(f(x,y)も同じ)を考え、その数値解を求めるプログラムを作成せよ。また、オイラー法と2次ルンゲクッタ法の実行結果を示して、2つの近似精度を比較せよ。 3 1と同じ常微分方程式(f(x,y)も同じ)を考え、その数値解を4次ルンゲクッタ法を使って求めるプログラムを作成せよ。また、オイラー法、2次ルンゲクッタ法、4次ルンゲクッタ法の実行結果を示して、3つの近似精度を比較せよ。 以上の3つです。 休んでいた自分が悪いのですが、ネットで調べてもよくわからなくて… わかる方、よろしくおねがいします…

  • オイラー法

    ①dx(t)/dt=1-x²(t),Δt=1/10,x(0)=0であるときx(0.3)の値をオイラー法とエクセルを用いてグラフにせよ オイラー法を使いエクセルで上記の微分方程式を解きたいのですがエクセルの使い方がわからず悩んでいます 教えてください、お願いします

  • オイラー法 Excel

    dx(t)/dt=f(t,x(t)),Δt=1/2,x(0)=0であるとき、x(4)の値をオイラー法を用いてExcelで示せ。ただしf(t,x(t))={0.8t(x(t)<1),-0.6t(x(t)≧1)である こちらの問題をExcelで解ける方、教えてください!!どうかお願いします!!

  • ルンゲクッタ法(RK4)で斜方投射の問題を解く

    下のように自由落下する物体の速度の変化を求めるプログラムをC言語で作成しました。これを応用して斜方投射された物体のx、y座標の「位置」の変化を求めるプログラムをルンゲクッタ法(RK4)で作成したいです。x,yについて2つの式を立てないといけないと思うのですがどのような式を立ててプログラムを組めかよいかよく分かりません。どなたか式とできればプログラムを教えて下さい。 k:空気抵抗、m:物体の質量 g:重力加速度としています。 #include <stdio.h> double yd(double t,double y){ double k=0.1; double m=0.1; double g=9.8; double r=g-((k*y)/m); return r; } double runge(double t,double y,double h){ double k1,k2,k3,k4,r; double h2=h/2.0; k1=yd(t,y); k2=yd(t+h2,y+(h2*k1)); k3=yd(t+h2,y+(h2*k2)); k4=yd(t+h,y+(h*k3)); r=y+(h/6.0)*((k1+(2.0*k2)+(2.0*k3)+k4)); return r; } void main(){ double y=0.0; double h=0.001; double t=0.0; int i; for(i=0;i<100;i++){ t+=h; y=runge(t,y,h); printf("%f %f\n",t,y); } }

  • 単振動の微分方程式を刻みhについてルンゲクッタで求める。

    m*d^2x/dt^2=-kx x(0)=1 dx/dt=0 というのが与えられて二階微分だから一階微分にするために dx(t)/dt=v(t) dv(t)/dt=-k*x(t)/mという式を立てました。オイラー法ではできたのですが2次、4次のルンゲクッタだとできません。どなたか回答お願いします。

  • 微分方程式をオイラー法でときたい

    高階常微分方程式 y” =f(x、y) y(0)=yo y’(0)=y’o この式と初期値でオイラー法を使って解きたいのですが... オイラー法を二回使えばよい、一階の連立方程式に なおせばよい。という意味がいまいちつかめません。 教えていただけると助かります。お願いします。

  • オイラー法を利用した空気抵抗を加味した斜方投射の式

    オイラー法を利用した空気抵抗を加味した斜方投射の式について 空気抵抗を加味した斜方投射の式について教えて下さい F=kv^2で空気抵抗が与えられているオイラー法を用いて速度の式を立てる場合 v(n+1)=v(n)-a×⊿t で表せられると思うのですがこの式のaを a=(m×g-k×v(n)^2)/m で計算することは正しいのでしょうか 自分がこの式で計算したところ、yの値が減少せずに増加のみの状態になってしまいました また、終端速度がv=√(m×g/k)で表せられると思います 終端速度に達した場合vは等速直線運動になると思うのですがどのような式で表せばよいのかわかりません よろしくお願いします ※質量m=7.26、初速度v(0)=10、空気抵抗係数k=0.01、⊿t=0.1です 下に自分の考えたオイラー法のアルゴリズムを記述しておきます (u=x方向の加速度、v=y方向の加速度、、x=x方向の飛距離、y=y方向の飛距離) loop n=1 a=(m*g-k*v(n)^2)/m u(n+1)=u(n) v(n+1) = v(n) - a*⊿t x(n+1) = x(n) + u(n)*⊿t y(n+1) = y(n) + v(n)*⊿t end loop

  • ルンゲクッタ法で数値解析(C言語)

    コンピュータ(プログラミング)のカテゴリに投稿しようかとも考えましたが、物理学に関する数値計算ということなので、物理学のカテゴリに投稿しました。 単振動や減衰振動、そして強制振動などを数値解析したいと思い、ルンゲクッタ法を使いシミュレートしてみようとしています。ルンゲクッタ法という方法を全く知らなかったので、インターネットや図書館で調べたのですが、どうしても分からないことろがあるので質問することにしました。 書籍やネットの情報を参考にしながら、単振動の場合を数値解析してみました(C言語を使って)。この単振動はうまくできたのですが、どうしても、その先の、減衰振動の数値解析がうまくいかないので、困っています。 ---- #include <stdio.h> double f1(double t,double x,double v); double f2(double t,double x,double v); int main() { double t,x,v,dt,tmax; double k1[2],k2[2],k3[2],k4[2]; x=1.0; //位置の初期値 v=0.0; //速度の初期値 dt=0.01; //刻み幅 tmax=500.0; //繰り返し最大回数 FILE *output; output=fopen("output.dat","w"); for(t=0;t<tmax;t+=dt) { k1[0]=dt*f1(t,x,v); k1[1]=dt*f2(t,x,v); k2[0]=dt*f1(t+dt/2.0,x+k1[0]/2.0,v+k1[1]/2.0); k2[1]=dt*f2(t+dt/2.0,x+k1[0]/2.0,v+k1[1]/2.0); k3[0]=dt*f1(t+dt/2.0,x+k2[0]/2.0,v+k2[1]/2.0); k3[1]=dt*f2(t+dt/2.0,x+k2[0]/2.0,v+k2[1]/2.0); k4[0]=dt*f1(t+dt,x+k3[0],v+k3[1]); k4[1]=dt*f2(t+dt,x+k3[0],v+k3[1]); x=x+(k1[0]+2.0*k2[0]+2.0*k3[0]+k4[0])/6.0; v=v+(k1[1]+2.0*k2[1]+2.0*k3[1]+k4[1])/6.0; fprintf(output,"%f %f %f\n",t,x,v); } fclose(output); return 0; } double f1(double t,double x,double v) { return v; } double f2(double t,double x,double v) { return (-x); } ---- このソースは単振動のもので、減衰振動のときは、最後の方の ---- double f1(double t,double x,double v) { return v; } double f2(double t,double x,double v) { return (-x); ---- の部分が変わるのだと思うのですが、よく分かりません。 減衰振動は、mx"=-kx-rx'で表されるので、この式を変形(?)したものが、入るのかな、という予想程度にしか分かりません。 ソースをどのように変えれば減衰振動を解析できるのでしょうか。 どなたか詳しい方、教えてください。お願いします。