コードがなぜか発散してしまいます

Question

発生している問題・エラーメッセージ

t=0.14あたりで関数f1,f2が発散してしまう

解決したいこと

printfで発散する部分は見つけたのですがなぜ発散するのかがわかりません
関数f1とf2が発散してしまう理由を教えてくださるとありがたいです。

c
ソースコード

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

#define h 0.001
#define g 9.80665
#define pai 3.14159265359
#define R 2.04

double Euler(double y0,double y1,double y2,double y3,double t);
double f1(double c ,double x);
double f2(double z ,double c, double x);
double f3(double a);

int main()
{
Euler(1,0,2.04,1,0);

return 0;

}

double Euler(double y0,double y1,double y2,double y3,double t){
double x,z,a,c;
while(t<=5.0){
y0=f3(y0);
y1=f3(y1);
y2=f3(y2);
y3=f3(y3);
printf("%lf %lf %lf\n",Rcos(y0)sin(y2),Rsin(y0)sin(y2),Rcos(y0));
x=y0;
z=y1;
a=y2;
c=y3;
printf("%lf %lf %lf \n ",f1(c,x),f2(z,c,x),t);
y0=x+zh;
y2=a+c*h;
y1=z+f1(c,x)*h;
y3=c+f2(z,c,x)*h;
t=t+h;
//printf("%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf\n",x,z,a,y0,y4,y1,y5,t);

 }

}
double f1(double c ,double x){
return (Rccsin(x)cos(x)-gsin(x))/R;
}
double f2(double z ,double c, double x){
return -2zccos(x)/sin(x);
}
double f3(double a){
return a-pai;
}

Comments

[can110]

多分のコードブロック```～```の位置がおかしいためソースがコード表示されていないので修正ください。

[RaitoN]

質問が読みにくすぎて読む気にならない．

Answer 1

関数f1とf2が発散してしまう理由

t=0.14あたりでは、c=-2.73e+174、x=-4.10e+85 となっていますが、問題ありません?
f3()で、 a-pai としてますが、範囲を paiにしたいという事でしたら、引き算でなく、剰余(fmod) とすべきかと思います。
ただ、浮動小数点演算は必ず、誤差があるので注意を。

参考ですが、途中結果の出力、"%lf"より、"%e"とした方が、見やすいかと、、。

Comments

[linkinpark]

a-piは単振り子の運動なので角度をπからずれにしたいだけなので気にしないでください
申し訳ないのですがやはりf1が発散する理由はわからないのでしょうか？

[pepperleaf]

プログラムにコメントも無いので、意図することが分かりません。
y0, y1, y2, y3の値はそれぞれ、角度だと思うのですが、ループが回る度に大きく(正確には絶対値が大きく)なっています。sin(), cos()は、大きい値でも計算可能ですが、有効桁数の問題で誤差が大きくなります。±π の範囲にすべきでしょう。また、繰り返しの累積誤差も考慮されているとは思えません。
申し訳ないですが、これ以上は、きちんと誤差評価しないとちょっと分かりません。(ほぼ、計算誤差の問題と考えてます)
ちなみに、paiの桁数を少し増やした時の t=0.14での y0,y1,y2,y3は以下のようになりました。(手元の環境)
y0=-890.005, y1=-316532807.883, y2=1089.465, y3=-177204037.548
y0=-890.005, y1=-316532931.730, y2=1089.465, y3=-177204110.986

[linkinpark]

球面でのふりこの動きをオイラー法で解析したかったので誤差が出るのは大丈夫です。fmodを使って範囲を制限するといいんですねありがとうございます!

Answer 2

回答ではありません・・・治す気がなさそうなので整形しておきました。

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

#define h 0.001
#define g 9.80665
#define pai 3.14159265359
#define R 2.04

double Euler(double y0, double y1, double y2, double y3, double t);
double f1(double c, double x);
double f2(double z, double c, double x);
double f3(double a);

int main( )
{
    Euler(1, 0, 2.04, 1, 0);

    return 0;
}

double Euler(double y0, double y1, double y2, double y3, double t)
{
    double x, z, a, c;
    while(t <= 5.0) {
        y0 = f3(y0);
        y1 = f3(y1);
        y2 = f3(y2);
        y3 = f3(y3);
        printf("%lf %lf %lf\n", R * cos(y0) * sin(y2), R * sin(y0) * sin(y2), R * cos(y0));
        x = y0;
        z = y1;
        a = y2;
        c = y3;
        printf("%lf %lf %lf \n ", f1(c, x), f2(z, c, x), t);
        y0 = x + z * h;
        y2 = a + c * h;
        y1 = z + f1(c, x) * h;
        y3 = c + f2(z, c, x) * h;
        t  = t + h;
        //printf("%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf\n",x,z,a,y0,y4,y1,y5,t);
    }
}
double f1(double c, double x)
{
    return (R * c * c * sin(x) * cos(x) - g * sin(x)) / R;
}
double f2(double z, double c, double x)
{
    return -2 * z * c * cos(x) / sin(x);
}
double f3(double a)
{
    return a - pai;
}

double Euler()と言って置きながらreturnがないのは？・・・まぁ、main()でも受けてないから良いけど。

Answer 3

全然関係ないかもしれませんが…

例えば，位置エネルギーと運動エネルギーをやりとりするような話（バネとか）を
単純にオイラー法でシミュレートしようとすると，発散します．

バネの例で言えば，
位置を縦軸，速度を横軸にとった平面を考えましょう．
原点はエネルギーが０で，原点から遠ざかるほどエネルギーがでかい．
バネの(先端の)運動は，一定のエネルギー量を，位置←→運動　に行き来させていることになるので，
この平面上では円を描く形になります．

この円弧上での動きを，オイラー法というのは「円の接線方向にちょっと移動する」こととして近似する話であため，
その結果は円の外側にちょっと出てしまう．つまり，場のエネルギーがちょっと増えちゃう．
なので，何の工夫もしないと，イテレーション毎にエネルギーがどんどん増大していき，いつか発散します．

Comments

[linkinpark]

とてもわかりやすいたとえありがとうございます！オイラー法について理解が深まりました。ちなみに今回はπをもともとヘッダーにあったπにすると治りました！なぜかはまったくわかりません笑