気分転換(数値積分)

なんか今年度は試験がだいぶ早く、ここのところ試験対策を練りに練っています(寝ているということでは断じてありません)。

 

疲れた時にはアニメを観たり、自分が気になった事柄を調べたりしています。

 

要は気分転換ですね(^_^)

 

 

今日の気分転換は「数値積分」です。

 

数値積分をするうえで「区分求積法」というものが使われています。

 

区分求積法は数Ⅲ(?)で一応学習するのかな?覚えていませんが…(^_^;)

 

 

積分区間を細かく分割して長方形や台形に近似して積分していく方法です(ざっくりとした説明)。

 

例えば関数 f(x) を区間[a,b]で積分すると考える。このとき下図のように区間[a,b]をいくつかの区間に分けるんですね(下図は2等分しています)。

 

下図ではそれぞれの区間で長方形として面積を求めてそれらをすべて足していく、という感じです、はい。

 

このままでは積分値との誤差が大きいです。

 

これをもっと細かく分割していけば、積分値に近くなります。

f:id:Akyuine:20170711195015p:plain

 これをプログラミングで定積分を数値的に求めてみましょう。

 

矩形法、シンプソンの台形公式を使って定積分を求めたものと、今回はこれらとは少し異なる方法で数値積分してみようと思います。

 

それは分割した区間区間の両端と中点から最小二乗法を使って近似直線を求め、その区間での面積を出てきた近似直線を使って求めていくというものです。

ざっくり言えば下図のようにします。

f:id:Akyuine:20170711201957p:plain

今回は f(x) = x^3 を区間[0,1]で積分してみましょう。

 

これは積分値は手計算で求められますね(積分値は1/4)

 

これを先ほどの3つの方法で求めて、積分値との誤差の絶対値を表示するためのプログラムをC言語で作ってみました(区間を10分割しました)。

 

「integral_compare.c」

----------------------------------

#include <stdio.h>
#include <math.h>
#define N 10
#define Ans 0.25

//被積分関数
double f(double x)
{
return x*x*x;
}

//最小二乗法の傾き
double slope(double a,double b,double c)
{
double a1,a2;

a1=3*(a*f(a)+b*f(b)+c*f(c))-(a+b+c)*(f(a)+f(b)+f(c));
a2=3*(a*a+b*b+c*c)-(a+b+c)*(a+b+c);

return a1/a2;
}

//最小二乗法の切片
double segment(double a,double b,double c)
{
double b1,b2;

b1=(a*a+b*b+c*c)*(f(a)+f(b)+f(c))-(a*f(a)+b*f(b)+c*f(c))*(a+b+c);
b2=3*(a*a+b*b+c*c)-(a+b+c)*(a+b+c);

return b1/b2;
}

//各区間での面積
double minisize(double a,double b)
{
double g1,g2;

g1=slope(a,(a+b)/2,b)*a+segment(a,(a+b)/2,b);
g2=slope(a,(a+b)/2,b)*b+segment(a,(a+b)/2,b);

return (g1+g2)*(b-a)/2;
}

//矩形法
double integral1(double a,double b)
{
int k;
double h;
double s=0.0;

h=(b-a)/*1;
}
}

//それぞれの場合の理論値との誤差
int main(void)
{
//矩形法
printf("Pattern1: %.10f\n",fabs(Ans-integral1(0,1)));
//シンプソンの台形公式
printf("Pattern2: %.10f\n",fabs(Ans-integral2(0,1)));
//今回考えた求積法
printf("Pattern3: %.10f\n",fabs(Ans-integral3(0,1)));

return 0;
}

 

 -----------------------------------

実行結果

Pattern1: 0.0475000000
Pattern2: 0.0025000000
Pattern3: 0.0012500000

 

Pattern3は誤差が他の2つの方法に比べれば小さくなってますね

 

これだけで評価はできないので、分割数を増やして、分割数と積分値との誤差の関係を以下のグラフにまとめました。(gnuplot使用)

f:id:Akyuine:20170711204528p:plain

Pattern1が他の2つに比べて誤差が大きいですね。

 

なのでPattern2とPattern3だけを以下のグラフにまとめました。

f:id:Akyuine:20170711204542p:plain

気分転換でこんなに時間を使って、試験は大丈夫なのか自分で心配になってます…(^_^;)

 

まぁ今回考えた方法が自分が思っていたよりも積分値との誤差を小さくできていたので満足です。

 

*1:double)N);

for(k=0;k<N;k++)
{
s+=f(a+h*k)/N;
}

return s;
}

//シンプソンの台形公式
double integral2(double a,double b)
{
int k;
double h;
double s=0.0;

h=(b-a)/((double)N);

for(k=0;k<N;k++)
{
s+=(f(a+h*k)+f(a+h*(k+1)))/2/N;
}

return s;
}

//今回考えた求積法
double integral3(double a,double b)
{
int k;
double h;
double s=0.0;

h=(b-a)/((double)N);

for(k=0;k<N;k++)
{
s+=minisize(a+h*k,a+h*(k+1