２０２３年度前期 数値解析 第１１回授業ノート
２０２３年６月２９日 三浦憲二郎
スライド１ 講義アウトライン•

関数「」と「•「•「」」近似による関数近似」補間

スライド２ 関数近似 p.116•

複雑な関数を簡単な関数で近似する•

閉区間[a,b]で定義された関数 f(x)をｇ（ｘ）=∑aiφi(x)で近似関数「」する．関数系 φi(x) は[a, b]上で連続かつ１次独立関数が 1 次独立とは，関数の差の「」が最小になるように係数 ai を定める．
スライド３ 定理 6.1 p.117•(証明) (6.1)よりが最小になるには「」の「」条件を満たさなければならない．

連立１次方程式「」を解くことと同値．

1

スライド４ グラム行列式•

要素 x0, x1, …,xn-1 が線形独立であることと，以下の「（6.2）はただ「」行列式が 0 でないことと同値．

」の解を持つ．

定理 6.2 φi=xi (i=0,1,2,…)は [a,b] で連続であり，かつ 1 次独立である．
（証明） 任意の自然数 n に対して，と仮定すると，閉区間[a，ｂ]の任意の t に対して 0, したがって
なぜならば, n 次方程式の解は高々n 個．
スライド５ 関数近似の例

p.118

•

例 6.3 f(x)=x2 に対する最小 2 乗近似 1 次式 g(x)=a0 + a1 x を[0,1]で求めよ．

•

φ0(x)=1, φ1(x)=x とするとしたがって，a0 = -1/6, a1 = 1.

y=x-1/6.

よってスライド６ 定理 6.4 p.119
m 組の与えられたデータ(x0,y0), (x1,y1), … ,(xm-1,ym-1)を通る関数 f(x)を g(x)=∑ai φi(x)によって近似することを考える．
が最小になるように a0,a1,…,an-1 を決定する．

連立 1 次方程式の解である．

2スライド７ 最小 2 乗近似：プログラム 6.1 改p.120#include <stdio.h>#include <stdlib.h>#include <math.h>#define M 6/* データのペア数 */#define N 3/* N 次式で近似 *//* ベクトルの入力 */void input_vector2( double *b, char c, int n, FILE *fin, FILE *fout);/* 部分ピボット選択付きガウス消去法 */void gauss2( double a[N+1][N+1], double b[N+1], int n );/* 最小 2 乗近似 */void least_square( double *x, double *y, FILE *fout );int main(void){
FILE *fin, *fout;
double x[M], y[M];/* ファイルのオープン */if ( (fin = fopen( "input_func.dat", "r")) == NULL ){printf("ファイルが見つかりません : input_func.dat ¥n");exit(1);}if( (fout = fopen( "output_func.dat", "w")) == NULL ){printf("ファイルが作成できません : output_func.dat ¥n");exit(1);}
input_vector2( x, 'x', M, fin, fout );/* ベクトル x の入出力 */input_vector2( y, 'y', M, fin, fout );/* ベクトル y の入出力 */
least_square( x, y, fout ); /* 最小 2 乗近似 */fclose(fin); fclose(fout); /* ファイルのクローズ */return 0;}
void least_square( double x[M], double y[M], FILE *fout ){double a[N+1], p[N+1][N+1];
int i, j, k;/* 右辺ベクトルの作成 */for(i=0; i <= N; i++) {input_func.dat
a[i]=0.0;
0.0
0.2
0.4 0.6 0.8
1.0for( j = 0; j < M; j++) {a[i] += y[j]*pow(x[j],(double)i) ;2.0
2.12 1.62 2.57 1.53 2.0}}/* 係数行列の作成 */output_func.dat
for( i = 0; i <= N; i++) {ベクトル x は次の通りですfor( j = 0; j <= i; j++ ) {p[i][j]=0.0;
0.00
0.20
0.40
0.60
0.80
1.00
for( k =0; k < M; k++) {ベクトル y は次の通りです
p[i][j] += pow( x[k], (double)(i+j) );
}
2.00
2.12
1.62
2.57
1.53
2.00
p[j][i] = p[i][j];最小 2 乗近似式は}}
y = 0.38 x^3 - 0.76 x^2 + 0.28 x^1 + 2.00 x^0/* 連立 1 次方程式を解く. 結果は a に上書き */
gauss2( p, a, N+1 );/* 結果の出力 */3}fprintf( fout, "最小 2 乗近似式は y =");
for( i = N ; i >= 0 ; i--) {if(a[i]>0){
if(i==N){
fprintf(fout, " %5.2f x^%d ", a[i],i);}else{
fprintf(fout, "+ %5.2f x^%d ", a[i],i);
}
}
else{
fprintf(fout, "- %5.2f x^%d ", fabs(a[i]),i);}}
fprintf(fout, "¥n");

スライド８ ラグランジュ補間 p.124
点(x0,f0), (x1,f1), … , (xn,fn) が与えられたとき，これらすべての 点(xi,fi) を通る曲線 y=f(x) を求めてx0 < x < xn の与えられた点 以外の関数値を求めることを「」，「」 (interpolation)するとういう． fk = f(xk), k=0,1,…,n が与えられたとき，等式 P(xk) = fk, k=0,1,…,n を満たす多項式 P(x) をf(x) の「」式という．

定理 6.5 補間条件を満たす n 次多項式 Pn(x) はただ１つに定まる．

V の行列式：「
スライド９ ラグランジュ補間」(Vandermonde)の行列， 解が「p.125ｎ次のラグランジュ補間多項式，ラグランジュ補間(Lagrange interpolation)

基本多項式 li(x) の 点 xk (0≦k≦n) での値は4」に定まる．

スライド１０ ラグランジュ補間：プログラム
p.125#include <stdio.h>#include <stdlib.h>#define N 9/* ベクトルの入力 */void input_vector3( double b[N+1], char c, FILE *fin );/* ラグランジュ補間 */double lagrange( double x[N+1], double y[N+1], double xi );int main(void){
FILE *fin, *fout;
double x[N+1], y[N+1], xi; /* xiは補間点 */printf("補間点を入力してください--->");
scanf("%lf", &xi);/* ファイルのオープン */if ( (fin = fopen( "input_lag.dat", "r")) == NULL ){printf("ファイルが見つかりません : input_lag.dat ¥n");exit(1);}if( (fout = fopen( "output_lag.dat", "w")) == NULL ){printf("ファイルが作成できません : output_lag.dat ¥n");exit(1);}
input_vector3( x, 'x', fin );/* ベクトルxの入出力 */input_vector3( y, 'y', fin );/* ベクトルyの入出力 */printf("補間の結果は、P(%f)=%f¥n", xi, lagrange(x,y,xi) );}/* グラフを描くために結果をファイルに出力 */for( xi = x[0]; xi <= x[N]; xi += 0.01 ){fprintf(fout, "%f ¥t %f¥n", xi, lagrange(x,y,xi) );}fclose(fin); fclose(fout); /* ファイルのクローズ */return 0;

input_lag.dat
0.0

0.2

0.4

0.6

0.8

1.2

1.4

1.6

1.8

2.0

2.0

2.11.6

2.6

1.5

2.7

0.67

3.5

0.94

2.0

/* ラグランジュ補間 */double lagrange( double x[N+1], double y[N+1], double xi ){int i, k;double pn = 0.0, li;/* P_n(x)の計算 */for ( i = 0; i <= N ; i++ ){li = 1.0;/* l_i(x)の計算 */for( k = 0; k <= N; k++ ){if( k != i ) li *= (xi -x[k]) / (x[i]-x[k]);}pn += li * y[i];}}return pn;/* b[0...N]の入力 */void input_vector3( double b[N+1], char c, FILE *fin ){int i;
for( i = 0 ; i <= N ; i++){fscanf(fin, "%lf", &b[i]);}}5