ある振動の方程式をC言語を用いて、ルンゲクッタで解析するプログラムについての質...

[joecoal]

ある振動の方程式をC言語を用いて、ルンゲクッタで解析するプログラムについての質問です。

式で表すと、
x"+rx'+(k^2)sinx-((a^2)sin2(x-j))/4
となる振動モデルをルンゲクッタで解析します。
ｘはある一定の値に落ち着きます。

たとえばk=0.001、r=0.01、 A=0.005、 j=0 という条件でプログラムを実行すると、
xは徐々に増加し1.489に落ち着きます。（出力はラジアン）

しかし、ｊ に値を入れると出力ｘにマイナスの値が出力されます。
たとえば
ｊ＝10（度）のときは、ｘは-1.318...
ｊ＝50（度）のときは、ｘは-0.650...
が出力されてしまいます。
理想としては
ｊ＝10（度）のときは、ｘは0.190程度
ｊ＝50（度）のときは、ｘは0.567程度

具体的にどのように改良したらよいのでしょうか・・・＾＾；

プログラムの問題でしょうか？？
よろしくお願いします

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

#define NDIV 10000;

static double k=0.001;
static double r=0.01;
static double A=0.005;
static double j=30;

// dx/dt = kx
double f(double x,double v){
return v;
}
// dv/dt = kv
double g(double x,double v){//@
return -((r*v)+(k*k)*sin(x)-((A*A)*sin(2*(x-(j*(3.14159/180)))))/4);
}

int main(void){
double x; // 変位
double v; // 速度
double t; // 時刻
double h; // =dt
double hh; // =dt/2
double kx1,kv1;
double kx2,kv2;
double kx3,kv3;
double kx4,kv4;
double pi;
int i,n=NDIV;

FILE *fp1,*fp2,*fp3;

pi = acos(-1.);
h = 10;
hh = h/2;
x = 0.01;
v = 0;
fp1 = fopen("test-t.txt", "w");
fp2 = fopen("test-x.txt", "w");
fp3 = fopen("test-v.txt", "w");

if(fp1 == NULL){
printf("can't open\n");
}else{
printf("open\n");
}
printf("honban\n");
printf(" # t x v f(!=0)\n"); //A=0;
for(i=0;i<=n;i++){
t = h*i;
if(i%10==0){
printf("%3d %7.4f %9.4f %9.4f",i,t,x,v);
printf("\n");
}
fprintf(fp1, "%.7f\n",t);
fprintf(fp2, "%.7f\n",x);
fprintf(fp3, "%.7f\n",v);
kx1 = f(x,v);
kv1 = g(x,v);
kx2 = f(x+hh*kx1,v+hh*kv1);
kv2 = g(x+hh*kx1,v+hh*kv1);
kx3 = f(x+hh*kx2,v+hh*kv2);
kv3 = g(x+hh*kx2,v+hh*kv2);
kx4 = f(x+h*kx3,v+h*kv3);
kv4 = g(x+h*kx3,v+h*kv3);
x = x + h/6*(kx1+2*kx2+2*kx3+kx4);
v = v + h/6*(kv1+2*kv2+2*kv3+kv4);
}
fclose(fp1);
fclose(fp2);
fclose(fp3);
printf("honban3");
return 0;
}

[パコネコ]

私ではよくわかりませんが、ｆって必要ですか？
（その解析をよく知らないため）
ｆで渡してるｖをそのまま返すだけならｘに渡す必要ないですし・・・