２自由度減衰振動のルンゲクッタ法

２自由度の減衰振動をルンゲクッタ法で計算したいと思い、プログラムを組んでみました。
ただし、実行結果がどうも正しいものになっていないようです。
どなたか詳しい方、教えて下さらないでしょうか？
以下が組んでみたプログラムです。

#include <stdio.h>
#include <math.h>

#define NDIV 500;

static double m1=1.0;
static double m2=1.0;
static double k1=1.0;
static double k2=1.0;
static double c1=0.00;
static double c2=0.00;

// dx/dt = kx
double f(double x1,double v1){
return v1;
}
// dv/dt = kv
double g(double x1,double v1,double x2, double v2){
return -((c1+c2)*v1-c2*v2+(k1+k2)*x1-k2*x2)/m1;
}

double j(double x2, double v2){
return v2;
}

double p(double x1, double v1, double x2, double v2){
return (c2*v1-c2*v2+k2*x1-k2*x2)/m2;
}

int main(void){
double x1; // 変位
double x2;
double v1; // 速度
double v2;
double c1;
double c2; // 減衰
double t; // 時刻
double h; // =dt
double hh; // =dt/2
double kx11,kv11,kx21,kv21;
double kx12,kv12,kx22,kv22;
double kx13,kv13,kx23,kv23;
double kx14,kv14,kx24,kv24;
int i,n=NDIV;
char a[20],*pa;
FILE *fp1,*fp2,*fp3,*fp4,*fp5;


h = 0.1;
hh = h/2;
x1 = 1;
v1 = 0;
x2 = 0;
v2 = 0;
fp1 = fopen("ju-t.txt","w");
fp2 = fopen("ju-x1.txt","w");
fp3 = fopen("ju-v1.txt","w");
fp4 = fopen("ju-x2.txt","w");
fp5 = fopen("ju-v2.txt","w");

if(fp1 == NULL){
printf("can't open/n");
}
else{
printf("open/n");
}



printf(" # t x1 v1 x2 v2\n");
for(i=0;i<=n;i++){
t = h*i;
printf("i,t,x1,v1,x2,v2");
fprintf(fp1, "%.4f\n",t);
fprintf(fp2, "%.4f\n",x1);
fprintf(fp3, "%.4f\n",v1);
fprintf(fp4, "%.4f\n",x2);
fprintf(fp5, "%.4f\n",v2);
kx11 = f(x1,v1);
kv11 = g(x1,v1,x2,v2);
kx21 = j(x2,v2);
kv21 = p(x1,v1,x2,v2);
kx12 = f(x1+hh*kx11,v1+hh*kv11);
kv12 = g(x1+hh*kx11,v1+hh*kv11,x2+hh*kx21,v2+hh*kv21);
kx22 = j(x2+hh*kx21,v2+hh*kv21);
kv22 = p(x1+hh*kx11,v1+hh*kv11,x2+hh*kx21,v2+hh*kv21);
kx13 = f(x1+hh*kx12,v1+hh*kv12);
kv13 = g(x1+hh*kx12,v1+hh*kv12,x2+hh*kx22,v2+hh*kv22);
kx23 = j(x2+hh*kx22,v2+hh*kv22);
kv23 = p(x1+hh*kx12,v1+hh*kv12,x2+hh*kx22,v2+hh*kv22);
kx14 = f(x1+h*kx13,v1+h*kv13);
kv14 = g(x1+h*kx13,v1+h*kv13,x2+h*kx23,v2+h*kv23);
kx24 = j(x2+h*kx23,v2+h*kv23);
kv24 = p(x1+h*kx13,v1+h*kv13,x2+h*kx23,v2+h*kv23);
x1 = x1 + h*(kx11+2*kx12+2*kx13+kx14)/6;
v1 = v1 + h*(kv11+2*kv12+2*kv13+kv14)/6;
x2 = x2 + h*(kx21+2*kx22+2*kx23+kx24)/6;
v2 = v2 + h*(kv21+2*kv22+2*kv23+kv24)/6;
}
fclose(fp1);
fclose(fp2);
fclose(fp3);
fclose(fp4);
fclose(fp5);
return 0;
}

No.2 回答者： Tacosan 回答日時： 2014/06/19 10:34

「過去類似の研究を行っていた先輩が残してくれた実行結果と異なっていた」としても, パラメータが違っていたら当然違う形になりえるわけだからそれが直接「実行結果が間違っている」ことにはつながらないのでは?

で
http://oshiete.goo.ne.jp/qa/8642263.html
は放置ですか?

この回答へのお礼

パラメータは同じでしたが、先輩のとは前提条件が若干異なっていました。
今は、解決できました、回答ありがとうございました。

お礼日時：2014/06/20 16:54

No.1 回答者： Tacosan 回答日時： 2014/06/17 14:51

どういう微分方程式を解こうとしているのですか?

「実行結果がどうも正しいものになっていないようです」の根拠は?

この回答への補足

解こうとしている微分方程式は
m1*x1" = -(c1+c2)x1'+c2*x2'-(k1+k2)*x1+k2*x2
m2*x2" = c2*x1'-c2*x2'+k2*x1-k2*x2
の連立微分方程式です。
'は一階時間微分、"は二階時間微分のつもりです。
また、kはバネ係数、cは減衰係数です。

自分は大学生の研究室所属なのですが、過去類似の研究を行っていた先輩が
残してくれた実行結果と異なっていたので、正しくないのではと判断しました。

補足日時：2014/06/17 16:36