EMANの物理学 過去ログ No.7658 〜

 ● 連成振動のプログラミング

  投稿者:ヴィダル - 2009/09/14(Mon) 13:23  No.7658 
以前単振動のプログラミングについて質問させていただいたものです。
おかげさまで何とかうまくいきました!
次に多質点からなる系の代表格でもある連成振動についてプログラミングを書いているのですが・・・エラーばかりです(涙
どうやら配列がおかしいようですが・・・
調べたりしてもよく欠陥が分からないのでどうか教授ください><

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

#define dt 0.01
#define tmax 100
#define imax 128
#define m 1.00
#define k 1.00

/*アルゴリズム=4次のルンゲ=クッタ法=*/

int main(void){

double x,p,t;
double k1[2][imax],k2[2][imax],k3[2][imax],k4[2][imax];
int i;

/*初期値の設定*/

x[1]=1.00;
p[i]=0;

for(i=0;i<=imax;i++){

for(t=0;t<=tmax;t++){

k1[0]=dt*f1(t,x[i],p[i]);
k1[1]=dt*f2(t,x[i],p[i]);

k2[0]=dt*f1(t+dt/2.0,x[i]+k1[0]/2.0,p[i]+k1[1]/2.0);
k2[1]=dt*f2(t+dt/2.0,x[i]+k1[0]/2.0,p[i]+k1[1]/2.0);

k3[0]=dt*f1(t+dt/2.0,x[i]+k2[0]/2.0,p[i]+k2[1]/2.0);
k3[1]=dt*f2(t+dt/2.0,x[i]+k2[0]/2.0,p[i]+k2[1]/2.0);

k4[0]=dt*f1(t+dt,x[i]+k3[0],p[i]+k3[1]);
k4[1]=dt*f2(t+dt,x[i]+k3[0],p[i]+k3[1]);

x=x+(k1[0]+2.0*k2[0]+2.0*k3[0]+k4[0])/6;
p=p+(k1[1]+2.0*k2[1]+2.0*k3[1]+k4[1])/6;

}
printf("%f %f %f\n",t,x[i],p[i]);

}


/*関数の設定/*

double x,p;
int i;

double f1(double x, double p){

return(p[i]/m);

}

double f2(double x, double p){

return(-k*(2x[i]-x[i-i]-x[i+1]);

}

return 0;

}



  投稿者:全充 - 2009/09/14(Mon) 13:31  No.7659 
double x,p,t;

/*初期値の設定*/

x[1]=1.00;
p[i]=0;

xやp配列宣言してないけど、OKでしたっけ?

  投稿者:EMAN - 2009/09/14(Mon) 14:19  No.7660 
 関数の中で関数の定義があるってのも C 言語としては変ですね。
f1(double x, double p) と double f2(double x, double p) の定義は
main関数の外側へ出しちゃって下さい。

 そして、全充さんの言われる通り、
double x,p;
という宣言では配列は使えません。

double x[imax], p[imax];
などのようにしないと。

しかもこれらの変数と int i は、すべての関数で共通して使っているようですから、グローバル変数として関数の外、つまり main関数より前に
書いておかないといけません。

 それと、main関数から f1() とか f2() とかの関数を呼び出してますけど、
呼び出すときに引数が3つあるのに、定義されている関数では引数が二つしかありません。

 などなど、危険が一杯です。
 意味不明な箇所も一杯ありますので、
そういう部分は、具体的にこう直せ、と言いにくくなってます。

  投稿者:ヴィダル - 2009/09/14(Mon) 19:01  No.7664 
配列の宣言してなかった。。。ありがとうございます、直します。

>呼び出すときに引数が3つあるのに、定義されている関数では引数が二つしかありません。

ど、どういうことでしょうか?main関数の外に連立微分方程式の形を書けばOKということですか?

まずは、質点1(i=1を代入する)についてTMAX秒まで回して、それをIMAX回繰り返すという考えで作ったFOR文なのですが、この書き方じゃまずいですか?

  投稿者:全充 - 2009/09/14(Mon) 19:22  No.7665 
> k1[0]=dt*f1(t,x[i],p[i]);
> k1[1]=dt*f2(t,x[i],p[i]);
>
> k2[0]=dt*f1(t+dt/2.0,x[i]+k1[0]/2.0,p[i]+k1[1]/2.0);
> k2[1]=dt*f2(t+dt/2.0,x[i]+k1[0]/2.0,p[i]+k1[1]/2.0);
>
> k3[0]=dt*f1(t+dt/2.0,x[i]+k2[0]/2.0,p[i]+k2[1]/2.0);
> k3[1]=dt*f2(t+dt/2.0,x[i]+k2[0]/2.0,p[i]+k2[1]/2.0);
>
> k4[0]=dt*f1(t+dt,x[i]+k3[0],p[i]+k3[1]);
> k4[1]=dt*f2(t+dt,x[i]+k3[0],p[i]+k3[1]);

ここでf1やf2の呼び出しで引数が
f1(t,x,p)
と3つですが

関数宣言のほう

> double f1(double x, double p){
>
> return(p[i]/m);
>
> }
>
> double f2(double x, double p){
>
> return(-k*(2x[i]-x[i-i]-x[i+1]);
>
> }

f1(x,p)
になってます。
ということですね。

  投稿者:全充 - 2009/09/14(Mon) 19:37  No.7666 
よく見ると

> double f1(double x, double p){
>
> return(p[i]/m);
>
> }

m って何ですか?
x 使ってませんね

とか

> double f2(double x, double p){
>
> return(-k*(2x[i]-x[i-i]-x[i+1]);
>
> }

p 使ってませんよ
とか

とか

  投稿者:EMAN - 2009/09/15(Tue) 20:29  No.7669 
 昨晩はパソコンが使えない状況だったので、
昨日書いたのを書き込めませんでしたが、せっかく書いたので
投稿しておきます。

 元のプログラムに色々と不備があるので、
完璧に動く形にまで直せなかったし、
あまり改造すると元の形とは全く違ったものになるので、
とりあえず、こんな感じです。

http://homepage2.nifty.com/eman/etc/prog0915.txt

  投稿者:ヴィダル - 2009/09/15(Tue) 23:44  No.7670 
EMANさん、全充さんご指摘ありがとうございます。本当に丁寧で感謝×2です><

やっぱり一番は配列と変数の使い方がまずいってことですね。
もう一回コメントを参考にしながら頑張ってみます・・・!(今日も半日トライしていたのですが・・・w)