Simpson則 をテンプレートにして作成

開始行:
RIGHT:寄稿：東條遼平 * 定積分の近似 [#m0c9d72a] Simpson則(シンプソン則)を用いて定積分を求めてみます。 #ref(simpson1.png,nolink) このようなf(x)を積分区間[a,b]で積分するとします。 このとき十分大きな数nを考え、積分区間をn等分します。 PCでは無限を考えることができないために有限個にしか分割で... nを大きくすることで誤差を小さくすることができます。 ここでn等分された区間のうちの1つを拡大して見てみます。 #ref(simpson2.png,nolink) 分割された内の1つの区間[x0,x2]を考え、区間の幅をhとします。 [x0,x2]の中間値をx1とし、f(x0)、f(x1)、f(x2)をそれぞれ y0、y1、y2とします。 よく積分の説明で、矩形を考えて幅をゼロに近付けるというの... それができるのは無限という概念が存在するからで、有限の値... 扱えないPCで矩形を考えては大きなすき間が空いてしまい、 誤差が大きくなる、または分割数をかなり大きくしなければな... そのためこの(x0,y0)、(x1,y1)、(x2,y2)の三点を通る2次曲線... そうすることですき間を極めて小さくできると考えられます。 #ref(img87.png,nolink) とおいてx1を中心とした二次式を考えます。こうすると計算が... #ref(img88.png,nolink) #ref(img89.png,nolink) ここで、 #ref(img90.png,nolink) より、 #ref(img91.png,nolink) となり、 #ref(img92.png,nolink) これより[x0,x2]の範囲で積分すると、 #ref(img93.png,nolink) #ref(img94.png,nolink) #ref(img95.png,nolink) ここで積分区間をn分割したときのi番目の区間を積分した値は、 #ref(img96.png,nolink) #ref(img97.png,nolink) となり、積分区間[a,b]で積分すると、 #ref(img98.png,nolink) * ソースコード [#b891c478] 以上長くなりましたが、実装するのは簡単です。 double Simpson(double (*f)(double x), double x0, double ... { double h; double sum = 0.0, x; int i; h = (x1 - x0); if(num) h /= num; x = x0; for(i=0; i<num; i++){ sum += (*f)(x) + 4.*(*f)(x + h/2) + (*f)(x + h); x += h; } return sum*h/6; } 使いまわしがきくように関数ポインタを使用しました。 関数ポインタとは関数をアドレスで渡して使う方法で、動的に... という便利な方法です。main関数より、 double f1(double x) { return x; } double f2(double x) { return 4 / (1 + x * x); } という2つの関数を順に積分しています。関数名だけを書けば、... Simpson(f1, 2, 5, 1000); のようにすると、 for(i=0; i<num; i++){ sum += (*f)(x) + 4.*(*f)(x + h/2) + (*f)(x + h); x += h; } の部分は for(i=0; i<num; i++){ sum += f1(x) + 4.*f1(x + h/2) + f1(x + h); x += h; } と同じ処理をします。 - &ref(main.c); - &ref(calculation.c); - &ref(calculation.h);

終了行:
RIGHT:寄稿：東條遼平 * 定積分の近似 [#m0c9d72a] Simpson則(シンプソン則)を用いて定積分を求めてみます。 #ref(simpson1.png,nolink) このようなf(x)を積分区間[a,b]で積分するとします。 このとき十分大きな数nを考え、積分区間をn等分します。 PCでは無限を考えることができないために有限個にしか分割で... nを大きくすることで誤差を小さくすることができます。 ここでn等分された区間のうちの1つを拡大して見てみます。 #ref(simpson2.png,nolink) 分割された内の1つの区間[x0,x2]を考え、区間の幅をhとします。 [x0,x2]の中間値をx1とし、f(x0)、f(x1)、f(x2)をそれぞれ y0、y1、y2とします。 よく積分の説明で、矩形を考えて幅をゼロに近付けるというの... それができるのは無限という概念が存在するからで、有限の値... 扱えないPCで矩形を考えては大きなすき間が空いてしまい、 誤差が大きくなる、または分割数をかなり大きくしなければな... そのためこの(x0,y0)、(x1,y1)、(x2,y2)の三点を通る2次曲線... そうすることですき間を極めて小さくできると考えられます。 #ref(img87.png,nolink) とおいてx1を中心とした二次式を考えます。こうすると計算が... #ref(img88.png,nolink) #ref(img89.png,nolink) ここで、 #ref(img90.png,nolink) より、 #ref(img91.png,nolink) となり、 #ref(img92.png,nolink) これより[x0,x2]の範囲で積分すると、 #ref(img93.png,nolink) #ref(img94.png,nolink) #ref(img95.png,nolink) ここで積分区間をn分割したときのi番目の区間を積分した値は、 #ref(img96.png,nolink) #ref(img97.png,nolink) となり、積分区間[a,b]で積分すると、 #ref(img98.png,nolink) * ソースコード [#b891c478] 以上長くなりましたが、実装するのは簡単です。 double Simpson(double (*f)(double x), double x0, double ... { double h; double sum = 0.0, x; int i; h = (x1 - x0); if(num) h /= num; x = x0; for(i=0; i<num; i++){ sum += (*f)(x) + 4.*(*f)(x + h/2) + (*f)(x + h); x += h; } return sum*h/6; } 使いまわしがきくように関数ポインタを使用しました。 関数ポインタとは関数をアドレスで渡して使う方法で、動的に... という便利な方法です。main関数より、 double f1(double x) { return x; } double f2(double x) { return 4 / (1 + x * x); } という2つの関数を順に積分しています。関数名だけを書けば、... Simpson(f1, 2, 5, 1000); のようにすると、 for(i=0; i<num; i++){ sum += (*f)(x) + 4.*(*f)(x + h/2) + (*f)(x + h); x += h; } の部分は for(i=0; i<num; i++){ sum += f1(x) + 4.*f1(x + h/2) + f1(x + h); x += h; } と同じ処理をします。 - &ref(main.c); - &ref(calculation.c); - &ref(calculation.h);

ページ名: