２自由度減衰振動のルンゲクッタ法
２自由度の減衰振動をルンゲクッタ法で計算したいと思い、プログラムを組んでみました。
ただし、実行結果がどうも正しいものになっていないようです。
どなたか詳しい方、教えて下さらないでしょうか？
以下が組んでみたプログラムです。
#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;
}
お礼
ωを１として計算してみます！無事（？）計算出来ました。 説明有難う御座います。