Heun法 をテンプレートにして作成

開始行:
RIGHT:寄稿：東條遼平 * 常微分方程式の2次近似 [#aa49918b] 前回Euler法で Euler法による常微分方程式を解きました。 このとき得られる数値解において誤差が問題となりましたので もう少し正確に数値解を求めることが出来るHeun法を使ってみ... 実はこれでもあまり正確ではなく、一般的には4次近似のRunge-... Runge-Kutta法は次回作成してみたいと思います。 * Heun法とは [#v03eb99a] Euler法では次の図のように考えて数値解を逐次的に算出しまし... #ref(euler1.png,nolink) y(dt)を求める際にy'(0)のみを考慮していました。 これをHeun法ではy(dt)を求めるときの傾きをy'(0)だけではな... しています。とはいえ、y(dt)がわからなければy'(dt)もわかり... y'(0)を使って出てきたy(dt)を使ってy'(dt)を求め、y'(0)とy'... 傾きを計算し、改めてy(dt)を求めます。少しわかりにくいので... #ref(heun1.png,nolink) これはEuler法でy(2dt)まで求めたときのプロセスですが、この... y'(dt)とy'(0)の平均を使ってy(dt)を求めます。つまり次のよ... #ref(heun2.png,nolink) これは式で書くと、 #ref(img113.png,nolink) #ref(img114.png,nolink) #ref(img115.png,nolink) となりますが、k2nでfに渡している2つ目の引数を見てください。 この部分は上の図でいうとy(dt)に相当します。 また、k1はy'(0)に、k2はy'(dt)に相当するのがわかります。 直感的にも誤差が小さくなりそうです。 * ソースコード [#be3cdbf9] 前回Euler法を作っているためそれに少々手を加えるだけで内容... void Heun(double (*f[])(double t, double *x), double t0, double *x, double tn, int div, int ... { double h, t; double *k1, *k2, *temp; int i, j; if((k1 = (double*)malloc(sizeof(double)*num)) == NULL){ fprintf(stderr, "Allocation Error\n"); exit(1); } if((k2 = (double*)malloc(sizeof(double)*num)) == NULL){ fprintf(stderr, "Allocation Error\n"); exit(1); } if((temp = (double*)malloc(sizeof(double)*num)) == NUL... fprintf(stderr, "Allocation Error\n"); exit(1); } h = (tn - t0); if(div) h /= div; t = t0; for(i=0; i<div; i++){ for(j=0; j<num; j++){ k1[j] = (*f[j])(t, x); temp[j] = x[j] + h*k1[j]; } for(j=0; j<num; j++){ k2[j] = (*f[j])(t+h, temp); x[j] += (k1[j] + k2[j])*h/2; } t += h; } free(k1); free(k2); free(temp); } 一見ずいぶん長くなってるように見えますが、mallocとそのエ... 殆どです。 変数の数や次数を可変にしたために取り得る初期値の数が 一意に決まらない為、動的に領域を確保する必要がでてきます。 領域を確保してから引数で渡してもよかったのですが、そうす... 引数の数がやたらと増えてしまいますのでこのようにしました。 この方法だと、あるtにおけるy(t)の値だけを知りたいときには... グラフを作成するような場合(Euler法の記事の最後のソース)に... 延々と繰り返されますので無駄が大きいです。 - &ref(main.c); - &ref(calculation.c); - &ref(calculation.h);

終了行:
RIGHT:寄稿：東條遼平 * 常微分方程式の2次近似 [#aa49918b] 前回Euler法で Euler法による常微分方程式を解きました。 このとき得られる数値解において誤差が問題となりましたので もう少し正確に数値解を求めることが出来るHeun法を使ってみ... 実はこれでもあまり正確ではなく、一般的には4次近似のRunge-... Runge-Kutta法は次回作成してみたいと思います。 * Heun法とは [#v03eb99a] Euler法では次の図のように考えて数値解を逐次的に算出しまし... #ref(euler1.png,nolink) y(dt)を求める際にy'(0)のみを考慮していました。 これをHeun法ではy(dt)を求めるときの傾きをy'(0)だけではな... しています。とはいえ、y(dt)がわからなければy'(dt)もわかり... y'(0)を使って出てきたy(dt)を使ってy'(dt)を求め、y'(0)とy'... 傾きを計算し、改めてy(dt)を求めます。少しわかりにくいので... #ref(heun1.png,nolink) これはEuler法でy(2dt)まで求めたときのプロセスですが、この... y'(dt)とy'(0)の平均を使ってy(dt)を求めます。つまり次のよ... #ref(heun2.png,nolink) これは式で書くと、 #ref(img113.png,nolink) #ref(img114.png,nolink) #ref(img115.png,nolink) となりますが、k2nでfに渡している2つ目の引数を見てください。 この部分は上の図でいうとy(dt)に相当します。 また、k1はy'(0)に、k2はy'(dt)に相当するのがわかります。 直感的にも誤差が小さくなりそうです。 * ソースコード [#be3cdbf9] 前回Euler法を作っているためそれに少々手を加えるだけで内容... void Heun(double (*f[])(double t, double *x), double t0, double *x, double tn, int div, int ... { double h, t; double *k1, *k2, *temp; int i, j; if((k1 = (double*)malloc(sizeof(double)*num)) == NULL){ fprintf(stderr, "Allocation Error\n"); exit(1); } if((k2 = (double*)malloc(sizeof(double)*num)) == NULL){ fprintf(stderr, "Allocation Error\n"); exit(1); } if((temp = (double*)malloc(sizeof(double)*num)) == NUL... fprintf(stderr, "Allocation Error\n"); exit(1); } h = (tn - t0); if(div) h /= div; t = t0; for(i=0; i<div; i++){ for(j=0; j<num; j++){ k1[j] = (*f[j])(t, x); temp[j] = x[j] + h*k1[j]; } for(j=0; j<num; j++){ k2[j] = (*f[j])(t+h, temp); x[j] += (k1[j] + k2[j])*h/2; } t += h; } free(k1); free(k2); free(temp); } 一見ずいぶん長くなってるように見えますが、mallocとそのエ... 殆どです。 変数の数や次数を可変にしたために取り得る初期値の数が 一意に決まらない為、動的に領域を確保する必要がでてきます。 領域を確保してから引数で渡してもよかったのですが、そうす... 引数の数がやたらと増えてしまいますのでこのようにしました。 この方法だと、あるtにおけるy(t)の値だけを知りたいときには... グラフを作成するような場合(Euler法の記事の最後のソース)に... 延々と繰り返されますので無駄が大きいです。 - &ref(main.c); - &ref(calculation.c); - &ref(calculation.h);

ページ名: