プログラミング(c言語)

質問

以下のプログラムでオイラー法を用いた2階微分方程式を解きたいのですが、xとvyの値が変化しません。どこが間違っているのか教えて欲しいです
写真は問題文です

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

int main(){
int i;
double x,x0,y,y0,vx,vx0,vy,vy0,r,dt,t=0.0;
double G=6.674/pow(10.0,11.0),M=1.988*pow(10.0,30.0);
x0=1.496*pow(10.0,11.0),y0=0.0,vx0=0.0,vy0=2.978*pow(10.0,4.0);
dt=M_PI/100;
for(i=1;i<=2000;i++){
x=x0+dt*vx0; y=y0+dt*vy0;
r=sqrt(x0*x0+y0*y0);
vx=vx0-dt*(G*M*x0/(r*r*r));vy=vy0-dt*(G*M*y0/(r*r*r));
t=dt*i;
printf("x(%f)=%e y(%f)=%e vx(%f)=%e vy(%f)=%e\n",t,x,t,y,t,vx,t,vy);
x0=x; y0=y; vx0=vx; vy0=vy;
}
return(0);
}

回答No.1

asciiz · ベストアンサー

Gが10の-11乗オーダー
Mが10の30乗オーダー
x0は10の11乗オーダー
vy0が10の4乗オーダー
dtは10の-2乗オーダー

えーっと。

double は有効数字15桁ほどなので、15乗オーダー以上違う数値同士の演算が正確に出来ません。

1.00000000000001 に
0.00000000000001 を足すと
1.00000000000002 になりますが、

1.00000000000001 に
0.000000000000001 を足しても
1.00000000000001 のまま変化しない、と言うことがあったりします。
1.000000000000011 が、表現できないので。
(16桁目は不正確ながら行けたかも…でも17・18桁目となるともうランダムに?)

数式的には正しくても、C言語の演算精度(15桁)のせいで計算できていない、という気がします。
ご質問のプログラムを眺めてるだけなので、どこでそれが起こってるのか突き止めてませんが。
オーダーの近い所で計算を済ませられるよう、プログラミングに工夫が必要かもしれません。

お礼

なるほど、自分では全く思いつかないことでした。工夫して再チャレンジしてみます。