数値解析
2023年度前期 第5週 [2023年5月18日]

静岡大学
工学研究科機械工学専攻
計測情報講座創造科学技術大学院情報科学専攻
三浦 憲二郎講義アウトライン [5月18日]•数値積分•ニュートン・コーツ公式•台形公式•シンプソン公式•重積分

数値積分の必要性 p.135•初等関数（しょとうかんすう）とは、複素数を変数とする多項式関数・指数関数・対数関数主値の四則演算・合成によって表示できる関数である．三角関数や双曲線関数、そして両者の逆関数主値も初等関数と考えられる．
•微分•初等関数の導関数は必ず初等関数になる•積分•初等関数の積分は初等関数であらわされるとは限らない
例 フレネル積分

数値的に積分する以外方法がない！

ニュートン・コーツ公式 p.135•定積分を求めるために，分点a=x0<x1<x2<…<xn=bをとり，fk=f(xk)を通るラグランジュ補間多項式 Pn(x) を考える．ただし，それぞれの分点は等間隔h=(b-a)/nで並んでいるとする．Pn(x）は多項式なので以下のInを計算することができる．

この近似積分公式をn+1点のニュートン・コーツ公式(Newton-Cotes)と呼ぶ．

台形公式 p.135•ニュートン・コーツ公式 n=1 とする．h=b-a, x0=a, x1=b, x=a+sh

これを台形公式(Trapezoidal rule)という．

台形公式 p.136•実際の計算 区間[a,b]を n 等分．その分点 xｋ, 各小区間[xk,xk+1]に適用

これを複合台形公式，あるいは台形公式という．

台形公式：プログラム p.137#include <stdio.h>/* 関数の定義 */double func1(double x);
double func2(double x);/* 台形公式 */double trapezoidal( double a, double b, int n, double (*f)(double) );/* 関数の定義 */double func1(double x){
return( 2.0/(x*x) );}int main(void){
int n=100;

double func2(double x){
return( 4.0 / (1.0+x*x) );}printf("2.0/(x*x)を[1,2]で積分します。分割数は%dです¥n", n);printf("結果は%20.15fです¥n",trapezoidal(1.0, 2.0, n, func1) );}printf("4.0/(1+x*x)を[0,1]で積分します。分割数は%dです¥n", n);printf("結果は%20.15fです¥n",trapezoidal(0.0, 1.0, n, func2) );実行結果

return 0;

2.0/(x*x)を[1,2]で積分します。分割数は100です
結果は/* 台形公式 */double trapezoidal( double a, double b, int n, double (*f)(double) ){double T, h;
int i;
h = ( b - a ) /n ; /* 刻み幅の指定 *//* 台形公式 */T = ( (*f)(a) + (*f)(b) ) / 2.0;
for ( i = 1; i < n; i++) T += (*f)( a + i*h );T *= h;}return T;

1.000029166020881

4.0/(1+x*x)を[0,1]で積分します。分割数は100です。
結果は3.141575986923129

シンプソン公式 p.139•ニュートン・コーツ公式 n=2 とする．
h=(b-a)/2, x0=a, x1=a+h, x2=a+2h=bx=a+sh

これをシンプソン公式(Ｓｉｍｐｓｏｎ’s rule)という．

シンプソン公式 p.140•実際の計算 区間[a,b]を 2n 等分．その分点 xｋ, 各小区間[x2k,x2k+2]に適用

これを複合シンプソン公式，
あるいはシンプソン公式という．

シンプソン公式：プログラム p.140#include <stdio.h>/* 関数の定義 */double func1(double x);
double func2(double x);/* シンプソン公式 */double simpson( double a, double b, int n, double (*f)(double) );/* 関数の定義 */double func1(double x){
return( 2.0/(x*x) );}int main(void){
int n=50;
printf("2.0/(x*x)を[1,2]で積分します。分割数は%dです¥n", 2*n);printf("結果は%20.15fです¥n",simpson(1.0, 2.0, n, func1) );double func2(double x){
return( 4.0 / (1.0+x*x) );}printf("4.0/(1+x*x)を[0,1]で積分します。分割数は%dです¥n", 2*n);printf("結果は%20.15fです¥n",simpson(0.0, 1.0, n, func2) );

実行結果

}

return 0;/* シンプソン公式 */double simpson( double a, double b, int n, double (*f)(double) ){double S, h;
int i;
h = ( b - a ) /( 2.0*n ) ; /* 刻み幅の指定 *//* シンプソン公式 */S = ( (*f)(a) + (*f)(b) ) ;for ( i = 1; i < n; i++) {S += 4.0*(*f)( a + (2.0*i-1.0)*h ) + 2.0*(*f)( a + 2.0*i*h );}S += 4.0*(*f)( a + (2.0*n-1.0)*h );S *= h/3.0;}

return S;

2.0/(x*x)を[1,2]で積分します。分割数は100です
結果は1.000000002582389

4.0/(1+x*x)を[0,1]で積分します。分割数は100です。
結果は3.141592653589754

重積分 p.144
•重積分

•台形公式
•シンプソン公式

重積分 p.144•台形公式 [φ(xi),Ψ(xi)] ,φ(xi)=y0<y1<…<ym=Ψ(xi) とm等分,k=(ym-y0)/m•シンプソン公式[φ(xi),Ψ(xi)] ,φ(xi)=y0<y1<…<y2m=Ψ(xi) と2m等分, k=(y2m-y0)/2m重積分：プログラム p.145#include <stdio.h>#include <stdlib.h>/* 被積分関数の定義 */double func(double x, double y);/* yの積分区間 */double phi(double x);
double psi(double x);/* ベクトル領域の確保 */double *dvector(int i, int j);/* ベクトル領域の解放 */void free_dvector(double *a, int i);/* 重積分用の台形公式 */double trapezoidal2( double a, double b, int m, int n, double (*p)(double),double (*q)(double), double (*f)(double,double) );int main(void){
int n=20, m=20;
printf("8x^2+4y を x=[1,2], y=[2-x,x^2] で積分します ¥n");printf("xの分割数は%d, yの分割数%d, 結果は%15.10f¥n", m, n,
trapezoidal2( 1.0, 2.0, m, n, phi, psi, func ) );return 0;}/* 重積分用の台形公式 */double trapezoidal2( double a, double b, int m, int n, double (*p)(double),double (*q)(double), double (*f)(double,double) ){double T, h, k, *F, x, y1, y2;
int i, j;
F = dvector( 0, n );h = ( b - a ) /n ;/* 刻み幅の指定(x方向) *//* F_i の計算 */for ( i = 0; i <= n; i++ ) {x = a + i*h;y1 = (*p)(x); y2 = (*q)(x);k = ( y2 - y1 )/m;/* 刻み幅の指定(y方向) */F[i] = ( (*f)(x, y1) + (*f)(x, y2) ) / 2.0;
for ( j = 1; j < m; j++ ) F[i] += (*f)(x, y1+j*k);
F[i] *= k;}/* 積分の計算 */T = ( F[0] + F[n] ) / 2.0;
for ( i = 1; i < n; i++) T += F[i];T *= h;free_dvector( F, 0 );return T;}/* 被積分関数の定義 */double func(double x, double y) {return( 8.0*x*x + 4.0*y );}/* yの積分区間 */double phi(double x) {return( 2.0-x );}double psi(double x){
return( x*x );}double *dvector(int i, int j) /* a[i]～a[i+j]の領域を確保 */ {double *a;if ( (a=(double *)malloc( ((j-i+1)*sizeof(double))) ) == NULL ) {printf("メモリが確保できません(from dvector)¥n");exit(1);}
return(a-i);}
void free_dvector(double *a, int i){free( (void *)(a + i) ); /* (void *)型へのキャストが必要 */
}-

まとめ
•数値積分
•ニュートン・コーツ公式•台形公式•シンプソン公式•重積分