【掲示板利用宣言】

　次のフォームをすべてチェックしてからご利用ください。

　私は

　題名と投稿者名は具体的に書きます。
　課題の丸投げはしません。
　ソースの添付は「HTML変換ツール」で字下げします。
　返信の引用は最小限にします。
　環境（OSとコンパイラ）や症状は具体的に詳しく書きます。
　返信の付いた投稿は削除しません。
　マルチポスト（多重投稿）はしません。

掲示板２

管理者用メニュー    ツリーに戻る    携帯用URL    ホームページ    ログ    タグ一覧
No.23801
ルンゲクッタ法のプログラム
投稿者---ハチミツ(2005/10/23 16:14:48)


d^2y/dx^2=-0.25dy/dx-1.25y+0.5という微分方程式をルンゲクッタ法を使って解を求めるというプログラムを作ったのですが、うまく動きません
どこが間違っているのか分かりません。間違っている所を教えてくださいお願いします。

以下に作成したプログラムを示します。

#include <stdio.h>

double dy(double x,double y,double z){
return z; /* dy/dx = z */
}
double dz(double x,double y,double z){
return -1.25*y-0.25*z+0.5; /* dz/dx = -1.25y-0.25z+0.5 */
}
int main(void)
{
double x,y,z;
double k[4],l[4];
double dt = 2.0;
double xs = 0.0; /* xの初期値の定義 */
double ys = 0.0; /* yの初期値の定義 */
double zs = 0.1; /* y'の初期値の定義 */
y = ys;
z = zs;

for(x=0.0;x<=40.0;x+=0.0) {

k[0] = dy(x,y,z); /* kの値を求める */
l[0] = dz(x,y,z);
k[1] = dy(x+dt/2,y+k[0]*dt/2,z+l[0]*dt/2);
l[1] = dz(x+dt/2,y+k[0]*dt/2,z+l[0]*dt/2);
k[2] = dy(x+dt/2,y+k[1]*dt/2,z+l[1]*dt/2);
l[2] = dz(x+dt/2,y+k[1]*dt/2,z+l[1]*dt/2);
k[3] = dy(x+dt/2,y+k[2]*dt/2,z+l[2]*dt/2);
l[3] = dz(x+dt/2,y+k[2]*dt/2,z+l[2]*dt/2);
x+=dt;
y=y+(dt/6)*(k[0]+2*k[1]+2*k[2]+k[3]);
z=z+(dt/6)*(l[0]+2*l[1]+2*l[2]+l[3]);

printf("%f %f \n",x,y);
}
return 0;
}


この投稿にコメントする

削除パスワード
発言に関する情報 題名 投稿番号 投稿者名 投稿日時
<子記事> Re:ルンゲクッタ法のプログラム 23810 Craft 2005/10/24 04:31:55
<子記事> Re:ルンゲクッタ法のプログラム 23814 たか 2005/10/24 19:07:06

No.23810
Re:ルンゲクッタ法のプログラム
投稿者---Craft(2005/10/24 04:31:55)


>d^2y/dx^2=-0.25dy/dx-1.25y+0.5という微分方程式をルンゲクッタ法を使って解を求めるというプログラムを作ったのですが、うまく動きません
>どこが間違っているのか分かりません。間違っている所を教えてくださいお願いします。

ルンゲクッタ法というのがなにかよくわかっていませんが、
http://www.wakayama-u.ac.jp/~miw/ip3/04/10/10.html
ここの連立常微分方程式〜という例題がヒントになりませんか？
非常によく似た式みたいなので、何が間違ってるかわかるかもしれません。

上記ページをみた限りでは、ルンゲクッタ法で解を求める場合、
非常に細かい定義があるみたいですが、掲示されたコードには、
各ステートメントが何をしているのかもかいてないので、
これだと誰も回答できないと思います。


この投稿にコメントする

削除パスワード
No.23814
Re:ルンゲクッタ法のプログラム
投稿者---たか(2005/10/24 19:07:06)


>どこが間違っているのか分かりません。間違っている所を教えてくださいお願いします。

刻み値dtが大きすぎると思います。それから、ループ変数に
xを使うのも変です。ループ変数は別に取るべきです。そうで
ないと、xに刻み値dtを加算した値が無視されてしまいます。

double dy(double x, double y, double z) {
  return z; /* dy/dx = z */
}
double dz(double x, double y, double z) {
  return -1.25 * y - 0.25 * z + 0.5; /* dz/dx = -1.25y-0.25z+0.5 */
}
int main(void)
{
  int i;
  double x, y, z;
  double k[4], l[4];
  double dt = 0.05;
  double xs = 0.0; /* xの初期値の定義 */
  double ys = 0.0; /* yの初期値の定義 */
  double zs = 0.1; /* y'の初期値の定義 */

  x = xs;
  y = ys;
  z = zs;

  for (i = 0; i <= 100; i++) {
    k[0] = dy(x, y, z); /* kの値を求める */
    l[0] = dz(x, y, z);
    k[1] = dy(x + dt / 2, y + k[0]*dt/2,z+l[0]*dt/2);
    l[1] = dz(x + dt / 2, y + k[0]*dt/2,z+l[0]*dt/2);
    k[2] = dy(x + dt / 2, y + k[1]*dt/2,z+l[1]*dt/2);
    l[2] = dz(x + dt / 2, y + k[1]*dt/2,z+l[1]*dt/2);
    k[3] = dy(x + dt / 2, y + k[2]*dt/2,z+l[2]*dt/2);
    l[3] = dz(x + dt / 2, y + k[2]*dt/2,z+l[2]*dt/2);
    x += dt;
    y += (dt / 6) * (k[0] + 2 * k[1] + 2 * k[2] + k[3]);
    z += (dt / 6) * (l[0] + 2 * l[1] + 2 * l[2] + l[3]);

    printf("%f %f \n", x, y);
  }
  return 0;
}



この投稿にコメントする

削除パスワード
No.23816
Re:ルンゲクッタ法のプログラム
投稿者---たか(2005/10/24 19:33:16)


自己レスですが、それか、xをループ変数として使うままに
するには、次のように刻み値dtを加算する部分を外に出して
しまうかですね。見通しが悪くなりますが。

#include <stdio.h>

double dy(double x, double y, double z) {
  return z; /* dy/dx = z */
}
double dz(double x, double y, double z) {
  return -1.25 * y - 0.25 * z + 0.5; /* dz/dx = -1.25y-0.25z+0.5 */
}
int main(void)
{
  double x, y, z;
  double k[4], l[4];
  double dt = 0.05;
  double xs = 0.0; /* xの初期値の定義 */
  double ys = 0.0; /* yの初期値の定義 */
  double zs = 0.1; /* y'の初期値の定義 */

  /* x = xs; */
  y = ys;
  z = zs;

  for (x = xs; x <= 40.0; x += dt) {
    k[0] = dy(x, y, z); /* kの値を求める */
    l[0] = dz(x, y, z);
    k[1] = dy(x + dt / 2, y + k[0] * dt / 2, z + l[0] * dt / 2);
    l[1] = dz(x + dt / 2, y + k[0] * dt / 2, z + l[0] * dt / 2);
    k[2] = dy(x + dt / 2, y + k[1] * dt / 2, z + l[1] * dt / 2);
    l[2] = dz(x + dt / 2, y + k[1] * dt / 2, z + l[1] * dt / 2);
    k[3] = dy(x + dt / 2, y + k[2] * dt / 2, z + l[2] * dt / 2);
    l[3] = dz(x + dt / 2, y + k[2] * dt / 2, z + l[2] * dt / 2);

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

    printf("%f %f \n", x + dt, y);
  }
  return 0;
}



この投稿にコメントする

削除パスワード
管理者用メニュー    ツリーに戻る    携帯用URL    ホームページ    ログ    タグ一覧