ルンゲクッタ法で数値解析（Ｃ言語）

締切済

質問者：risyu
質問日時：2008/07/11 02:08
回答数：2件

コンピュータ（プログラミング）のカテゴリに投稿しようかとも考えましたが、物理学に関する数値計算ということなので、物理学のカテゴリに投稿しました。

単振動や減衰振動、そして強制振動などを数値解析したいと思い、ルンゲクッタ法を使いシミュレートしてみようとしています。ルンゲクッタ法という方法を全く知らなかったので、インターネットや図書館で調べたのですが、どうしても分からないことろがあるので質問することにしました。

書籍やネットの情報を参考にしながら、単振動の場合を数値解析してみました（Ｃ言語を使って）。この単振動はうまくできたのですが、どうしても、その先の、減衰振動の数値解析がうまくいかないので、困っています。
----
#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'で表されるので、この式を変形（？）したものが、入るのかな、という予想程度にしか分かりません。
ソースをどのように変えれば減衰振動を解析できるのでしょうか。

どなたか詳しい方、教えてください。お願いします。

通報する

この質問への回答は締め切られました。

質問の本文を隠す

回答 (2件)

最新から表示
回答順に表示

No.2

回答者： Tacosan
回答日時：2008/07/11 18:07

そんだけわかっているんだから, そのまま書けばいいだけだろうに....

単振動の場合は mx'' = -kx で, v = dx/dt を導入することで
dx/dt = v, m dv/dt = -kx
という連立 1階微分方程式にしたんですよね? 減衰振動の mx'' = -kx - rx' も同じように連立 1階微分方程式にするだけです. つまり, こちらでも v = dx/dt を導入して
dx/dt = f1(t, x, v), dv/dt = f2(t, x, v)
の形にするだけ.