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

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

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;
}

8524ba23

2020/07/10 19:33 編集

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

RaitoN

2020/07/11 00:36

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

行動規範の内容に同意します

回答3件

ベストアンサー

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

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

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

投稿2020/07/11 02:35

pepperleaf

総合スコア6385

linkinpark

2020/07/11 04:17

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

pepperleaf

2020/07/11 06:36

プログラムにコメントも無いので、意図することが分かりません。 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

2020/07/11 07:35

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

行動規範の内容に同意します

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

c
1#include <stdio.h>
2#include <math.h>
3
4#define h 0.001
5#define g 9.80665
6#define pai 3.14159265359
7#define R 2.04
8
9double Euler(double y0, double y1, double y2, double y3, double t);
10double f1(double c, double x);
11double f2(double z, double c, double x);
12double f3(double a);
13
14int main( )
15{
16    Euler(1, 0, 2.04, 1, 0);
17
18    return 0;
19}
20
21double Euler(double y0, double y1, double y2, double y3, double t)
22{
23    double x, z, a, c;
24    while(t <= 5.0) {
25        y0 = f3(y0);
26        y1 = f3(y1);
27        y2 = f3(y2);
28        y3 = f3(y3);
29        printf("%lf %lf %lf\n", R * cos(y0) * sin(y2), R * sin(y0) * sin(y2), R * cos(y0));
30        x = y0;
31        z = y1;
32        a = y2;
33        c = y3;
34        printf("%lf %lf %lf \n ", f1(c, x), f2(z, c, x), t);
35        y0 = x + z * h;
36        y2 = a + c * h;
37        y1 = z + f1(c, x) * h;
38        y3 = c + f2(z, c, x) * h;
39        t  = t + h;
40        //printf("%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf\n",x,z,a,y0,y4,y1,y5,t);
41    }
42}
43double f1(double c, double x)
44{
45    return (R * c * c * sin(x) * cos(x) - g * sin(x)) / R;
46}
47double f2(double z, double c, double x)
48{
49    return -2 * z * c * cos(x) / sin(x);
50}
51double f3(double a)
52{
53    return a - pai;
54}
55