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

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

締切り済みの質問

n自由度の振動のルンゲクッタ法

n自由度の振動をルンゲクッタ法で計算したいと思い、プログラムの練習をしています。
現在2自由度まで計算するプログラムは組めているのですが、
nを入力し、そのnに対応した計算を行うルンゲクッタのプログラムが分かりません。
どなたか詳しい方に教えて頂きたいです。
以下に2自由度のプログラムを表示します。

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

#define NDIV 500;

static double m1=1.0;
static double m2=1.00;
static double k1=1.0;
static double k2=1.0;
static double c1=0.30;
static double c2=0.30;

// 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-23 20:26:19

QNo.8650584

すぐに回答ほしいです

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

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

回答(1)

ANo.1

それじゃあ,まず変数とかを書き並べるんじゃなくて配列を使って書き直してみようよ。

投稿日時 - 2014-06-23 22:01:14

お礼

そうですね、まずはそうしてみます。

投稿日時 - 2014-06-23 22:23:10

あなたにオススメの質問