C言語関係掲示板

過去ログ

No.968 ルンゲクッタの2階微分方程式のプログラム

[戻る] [ホームページ]
No.12477

ルンゲクッタの2階微分方程式のプログラムです。
投稿者---困ってます(2004/02/04 22:52:47)


yの値がうまく出せません。
#include <stdio.h>
#include <math.h>
#include <stdlib.h>


double dy(double x,double y,double z){ /* 導関数の定義 */
return z; /* dy/dx = z */
}

double dz(double x,double y,double z){ /* 導関数の定義 */
return -8*y-z; /* dz/dx = -8y-z */
}

int main()
{
printf("\n** d^2y/dx^2+dy/dx+8y=0 **\n");
int n=0;
double x,y,z,h;
double k[4],l[4];
double hKizami = 0.05; /* Δxの定義 */
double xSyoki = 0.0; /* xの初期値の定義 */
double ySyoki = 0.5; /* yの初期値の定義 */
double zSyoki = 1.0; /* z(y')の初期値の定義 */

printf("x 0=%5.2f y 0=%f \n",n,x,n,y);

FILE *fo;
fo=fopen("runge2.dat","w");
if (fo==NULL){
printf("Can't Open File!");
exit(1);
}


for(x=0.0;x<=5.0;x+=0.0)
{
n++;

k[0] = dy(x,y,z); /* kの値を求める */
l[0] = dz(x,y,z);

k[1] = dy(x+hKizami/2,y+k[0]*hKizami/2,z+l[0]*hKizami/2);
l[1] = dz(x+hKizami/2,y+k[0]*hKizami/2,z+l[0]*hKizami/2);

k[2] = dy(x+hKizami/2,y+k[1]*hKizami/2,z+l[1]*hKizami/2);
l[2] = dz(x+hKizami/2,y+k[1]*hKizami/2,z+l[1]*hKizami/2);

k[3] = dy(x+hKizami,y+k[2]*hKizami,z+l[2]*hKizami);
l[3] = dz(x+hKizami,y+k[2]*hKizami,z+l[2]*hKizami);

x+=hKizami;
y=y+(hKizami/6)*(k[0]+2*k[1]+2*k[2]+k[3]);
z=z+(hKizami/6)*(l[0]+2*l[1]+2*l[2]+l[3]);

printf("x%4d=%5f y%4d=%f \n",n,x,n,y);
fprintf(fo,"x%4d=%4.1f y%4d=%f \n",n,x,n,y);
}
fclose(fo);
return 0;
}



No.12491

Re:ルンゲクッタの2階微分方程式のプログラムです。
投稿者---たか(2004/02/05 14:57:09)


C言語は変数の定義はブロックの先頭でしか行えません。
それと、yとzの初期値が与えられていません。これを直せば動きました。

ここも参考にして下さい。オイラー法はきざみ値が
大きくなるとルンゲ・クッタ法よりも大きな影響(誤差が大きくなる)を
受けます。ルンゲ・クッタ法はシンプソン則を元に作られています。

#include <stdio.h>
#include <stdlib.h>

double dy(double x,double y,double z){ /* 導関数の定義 */
return z; /* dy/dx = z */
}

double dz(double x,double y,double z){ /* 導関数の定義 */
return -8*y-z; /* dz/dx = -8y-z */
}

int main(void)
{
  int n=0;
  double x,y,z,h;
  double k[4],l[4];
  double hKizami = 0.05; /* Δxの定義 */
  double xSyoki = 0.0; /* xの初期値の定義 */
  double ySyoki = 0.5; /* yの初期値の定義 */
  double zSyoki = 1.0; /* z(y')の初期値の定義 */
  FILE *fo;

  y = ySyoki; z = zSyoki;
  printf("\n** d^2y/dx^2+dy/dx+8y=0 **\n");
  printf("x 0=%5.2f y 0=%f \n",n,x,n,y);

  fo=fopen("runge2.dat","w");
  if (fo==NULL){
    printf("Can't Open File!");
    exit(1);
  }

  for(x=0.0;x<=5.0;x+=0.0) {
    n++;

    k[0] = dy(x,y,z); /* kの値を求める */
    l[0] = dz(x,y,z);

    k[1] = dy(x+hKizami/2,y+k[0]*hKizami/2,z+l[0]*hKizami/2);
    l[1] = dz(x+hKizami/2,y+k[0]*hKizami/2,z+l[0]*hKizami/2);

    k[2] = dy(x+hKizami/2,y+k[1]*hKizami/2,z+l[1]*hKizami/2);
    l[2] = dz(x+hKizami/2,y+k[1]*hKizami/2,z+l[1]*hKizami/2);

    k[3] = dy(x+hKizami,y+k[2]*hKizami,z+l[2]*hKizami);
    l[3] = dz(x+hKizami,y+k[2]*hKizami,z+l[2]*hKizami);

    x+=hKizami;
    y=y+(hKizami/6)*(k[0]+2*k[1]+2*k[2]+k[3]);
    z=z+(hKizami/6)*(l[0]+2*l[1]+2*l[2]+l[3]);

    printf("x%4d=%5f y%4d=%f \n",n,x,n,y);
    fprintf(fo,"x%4d=%5f y%4d=%f \n",n,x,n,y);
  }
  fclose(fo);
  return 0;
}


No.12544

Re:ルンゲクッタの2階微分方程式のプログラムです。
投稿者---困ってます。(2004/02/07 01:14:50)


たかさん
どうもありがとうございました。