ルンゲクッタ法で数値解析（Ｃ言語）
コンピュータ（プログラミング）のカテゴリに投稿しようかとも考えましたが、物理学に関する数値計算ということなので、物理学のカテゴリに投稿しました。
単振動や減衰振動、そして強制振動などを数値解析したいと思い、ルンゲクッタ法を使いシミュレートしてみようとしています。ルンゲクッタ法という方法を全く知らなかったので、インターネットや図書館で調べたのですが、どうしても分からないことろがあるので質問することにしました。
書籍やネットの情報を参考にしながら、単振動の場合を数値解析してみました（Ｃ言語を使って）。この単振動はうまくできたのですが、どうしても、その先の、減衰振動の数値解析がうまくいかないので、困っています。
----
#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'で表されるので、この式を変形（？）したものが、入るのかな、という予想程度にしか分かりません。
ソースをどのように変えれば減衰振動を解析できるのでしょうか。
どなたか詳しい方、教えてください。お願いします。
お礼
的確な回答、誠にありがとうございました。 文献を参考にプログラムを書いてみます。