Transcript Document

プログラミング論
数値積分
http://www.ns.kogakuin.ac.jp/~ct13140/Prog/
L-1
積分
区分求積法と台形公式のみ解説
多くの変数が登場する．しっかり整理すること
L-2
積分
• 不定積分は，「微分」の逆操作．すなわち，
f(x)を不定積分する/の不定積分を求め
るとは，「微分したらf(x)になる関数を探
す」こと．
• 定積分は，関数f(x)とx軸に挟まれる部
分の面積を求めること．
• 本講義では，定積分のみを扱う．
L-3
定積分
b
•  f ( x)dx とは？
a
x=b
f(x)
x=a
y=f(x)と，x軸と
x=aと，x=bで
囲まれた領域の面積を
求めること
x
L-4
区分求積法
• 細かい区間毎の長方形の面積の合計で近似．
b
全短冊の面積の合計と， f ( x)dx はほぼ等しい
a
f(xi) f(xi+1)
f(x)
y=f(x)
この短冊の面積は，
(xi+1－xi) f(xi)
a
xi xi+1
x
b
L-5
区分求積法
• a=x0<x1<…<xn-1<xn=bとして，
f (x0)，f (x1)，…，f (xn)が既知なら
n 1
a f ( x)dx   f ( xk )(xk 1  xk )
b
k 0
L-6
f(xi) f(xi+1)
左の値を使用
xi xi+1
x
f(xi) f(xi+1)
f(x)
• 当然，
長方形の高さは
左の値を使用し
ても，
右の値を使用し
ても良い．
f(x)
区分求積法
右側の値を使用
xi xi+1
x
L-7
区分求積法の誤差
f(xi+1)
f(x)
f(xi)
この部分が
近似の誤差
xi
xi+1
x
L-8
区分求積法
• より細かい区間に分割するほど誤差は少なくな
る．
• 長方形でなく台形で近似した方が，より正確な近
似になると期待される．
– 台形公式 (これは1次近似である)
• 近傍3点を用い，2次式で近似した方がより正確
な近似になると期待される．
– シンプソンの公式
• 3次式，4次式で近似した方が更に正確．
– ニュートン・コーツの積分公式
L-9
区分求積法台形公式
• 長方形ではなく，台形で近似する．
– より少ない誤差で計算が可能
f(xi)
f(x)
f(xi+1)
xi
x
xi+1
L-10
区分求積法台形公式
• a=x0<x1<…<xn-1<xn=bとして，
f (x0)，f (x1)，…，f (xn)が既知なら
n 1
1
a f ( x)dx   2  f ( xk )  f ( xk 1)( xk 1  xk )
k 0
b
ba
• 区間[a,b]をn等分した場合，h 
として
n
n

1
b
h

f ( x)dx 
f ( xk )  f ( xk 1 )
a
2 k 0


L-11
C言語プログラミング
b
•  f ( x)dx を求めるには？
a
h = (b-a)/N;
N当分すると仮定
f(x)
f(x1) f(x2)
h
a x 1 x2 x3 x4 x
x0
bx
n
L-12
C言語プログラミング
b
•  f ( x)dx を求めるには？
a
x0～xnを表示してみる．
N当分すると仮定
f(x1) f(x2)
f(x)
h = (b-a)/N;
x0 = a;
x1 = a + h*1;
x2 = a + h*2;
:
xi = a + h*i;
なので，
h
for(i=0; i<=N; i++){
printf("%lf\n", a+h*i); ax0 x1 x2 x3 x4
}
x
bx
n
L-13
C言語プログラミング
b
•  f ( x)dx を求めるには？
a
N当分すると仮定
各短冊の面積を表示してみる．
xi～xi+1の短冊の横の長さは h，
縦の長さは f (xi) なので，
h = (b-a)/N;
f(x)
for(i=0; i<N; i++){
xi = a+h*i;
f(x1) f(x2)
h
area = f(xi)*h;
printf("%lf\n",area);
}
関数 double f(double d)が
存在していると仮定.
ax x1 x2 x3 x4
0
x
bx
n
L-14
C言語プログラミング
b
•  f ( x)dx を求めるには？
a
全短冊の面積の合計を出す．
すなわち，定積分．
N当分すると仮定
f(x1) f(x2)
for(i=0; i<N; i++){
xi = a+h*i;
f(x)
h = (b-a)/N;
sum = 0.0;
h
area = f(xi)*h;
sum += area;
}
printf("%lf\n", sum);
ax x1 x2 x3 x4
0
x
bx
n
L-15
例(伊)
• y=x2 と x軸と x=1 と x=3 で囲まれ
た領域の面積は?
12
y=x2
10
8
y 6
4
2
0
0
0.5
1
1.5
x
2
2.5
3
3.5
L-16
例(伊)
100分割，非台形
#define NUM 100
#include <stdio.h>
void main(){
int i;
double x, y, sum, h;
h = 2.0/NUM;
sum = 0.0;
for(i=0; i<NUM; i++){
x = 1.0 + h*i;
y = x*x;
sum += y*h;
}
printf("%lf\n", sum);
}
L-17
#define NUM 1000
#include <stdio.h>
void main(){
int i;
double xa, xb, ya, yb, sum, h;
h = 2.0/NUM;
sum = 0.0;
for(i=0; i<NUM; i++){
xa = 1.0 + h*i;
xb = 1.0 + h*(i+1);
ya = xa*xa;
yb = xb*xb;
sum += (ya+yb)*h/2.0;
}
printf("%lf\n", sum);
}
例(伊)
1000分割，台形
L-18
例(伊)
3 2
26
x dx 
 8.66666 
1
3

100分割
1000分割
面積
誤差
面積
誤差
左値を
使用
8.5868
0.079867
8.658668
0.007999
右値を
使用
8.7468
0.080133
8.674668
0.008001
台形
8.6668
0.000133
8.666668
0.000001
L-19
例(呂)
• 下の扇形(1/4円)の面積は?
 x, 1  x 2 




1.0
0.8
0.6
0.4
0.2
0
0
0.2
0.4
0.6
0.8
1.0
L-20
例(呂)
1000分割，非台形
#define NUM 1000
#include <stdio.h>
#include <math.h>
double f(double x){
return
sqrt(1-x*x);
}
void main(){
int i;
面積=0.785889
double x, y, sum, h;
PI=3.143555
h = 1.0/NUM;
sum = 0.0;
for(i=0; i<NUM; i++){
x = h*i;
y = f(x);
sum += y*h;
}
printf("面積=%lf\nPI=%lf\n", sum, sum*4);
}
L-21