みんなの「教えて(疑問・質問)」にみんなで「答える」Q&Aコミュニティ

こんにちはゲストさん。会員登録(無料)して質問・回答してみよう!

締切り済みの質問

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

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

#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;
}

投稿日時 - 2014-06-17 12:46:20

QNo.8641922

すぐに回答ほしいです

このQ&Aは役に立ちましたか?

1人が「このQ&Aが役に立った」と投票しています

回答(2)

ANo.2

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


http://okwave.jp/qa/q8642263.html
は放置ですか?

投稿日時 - 2014-06-19 10:34:25

お礼

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

投稿日時 - 2014-06-20 16:54:29

ANo.1

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

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

投稿日時 - 2014-06-17 14:51:36

補足

解こうとしている微分方程式は
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:04

あなたにオススメの質問