三浦 憲二郎
静岡大学
創造科学技術大学院
情報科学専攻
工学部機械工学科
計測情報講座
数値解析
2019年度前期 第13週 [7月11日]
講義アウトライン [7月11日]
•関数近似と補間
•最小2乗近似による関数近似
•ラグランジュ補間
関数近似 p.116
•複雑な関数を簡単な関数で近似する
関数近似
•閉区間[a,b]で定義された関数f(x)をｇ（ｘ）=∑a
i
φ
i
(x)で近似
する．関数系 φ
i
(x) は[a, b]上で連続かつ１次独立
関数が1次独立とは，
関数の差のノルムの二乗
が最小になるように係数 a
i
を定める．
関数近似 確認問題
1. 2点(1,-1), (3,7)を通る直線を求めよ。
ヒント
y=ax+b
2. 次の関数の最小値を求めよ．
𝑓𝑓 𝑥𝑥, 𝑦𝑦 = 𝑥𝑥
2
+ 𝑥𝑥 𝑦𝑦 + 1 + 𝑦𝑦
2
3.最小2 乗法を用いて，3 点(1,-1), (2, 1), (3, 7) を近似する
直線を求めよ．
関数近似 確認問題 解答
1. 2点(1,-1), (3,7)を通る直線を求めよ。
(y+1)=a(x-1), 7+1=a(3-1), a=4, b=-5
2.次の関数の最小値を求めよ．
𝑑𝑑𝑑𝑑
𝑑𝑑𝑑𝑑
= 2𝑥𝑥 + 𝑦𝑦 + 1 = 0
,
𝑑𝑑𝑑𝑑
𝑑𝑑𝑦𝑦
= 𝑥𝑥 + 2𝑦𝑦 = 0
𝑥𝑥 =
2
3
, 𝑦𝑦 = −
1
3
,
3.最小2 乗法を用いて，3 点(1,-1), (2, 1), (3, 7) を近似する
直線を求めよ．
𝑦𝑦 = 4𝑎𝑎 −
17
3
定理6.1 p.117
•関数系 φ
i
が１次独立であれば，(6.1)のノルム||・||
2
2
を最小にする係数は
連立1次方程式
の解として一意に決定される．ここで，[a, b] 上の連続関数f(x)とg(x)に対して
(𝜑𝜑
0
, 𝜑𝜑
0
)
⋯
(𝜑𝜑
0
, 𝜑𝜑
𝑛𝑛−1
)
⋮
⋱
⋮
(𝜑𝜑
𝑛𝑛−1
, 𝜑𝜑
0
) ⋯ (𝜑𝜑
𝑛𝑛−1
, 𝜑𝜑
𝑛𝑛−1
)
𝑎𝑎
0
⋮
𝑎𝑎
𝑛𝑛−1
=
(𝑓𝑓, 𝜑𝜑
0
)
⋮
(𝑓𝑓, 𝜑𝜑
𝑛𝑛−1
)
定理6.1 p.117
•(証明) (6.1)より
が最小になるには極値の必要条件を満たさなければならない．
グラム行列式 p.117
•要素 x
0
, x
1
, …,x
n-1
が線形独立であることと，以下のグラム行列式が0でない
ことと同値．（6.2）はただ１つの解を持つ．
定理6.2 φ
i
=x
i
(i=0,1,2,…)は [a,b] で連続であり，かつ1次独立である．
（証明） 任意の自然数 n に対して，
と仮定すると，閉区間[a，ｂ]の任意の t に対して0, したがって
なぜならば, n次方程式の解は高々n個．
関数近似の例 p.118
•例6.3 f(x)=x
2
に対する最小2乗近似1次式g(x)=a
0
+ a
1
x を[0,1]で求めよ．
•φ
0
(x)=1, φ
1
(x)=x とすると
定理6.4 p.119
•m組の与えられたデータ(x
0
,y
0
), (x
1
,y
1
), … ,(x
m-1
,y
m-1
)を通る関数f(x)を
g(x)=∑a
i
φ
i
(x)によって近似することを考える．
が最小になるようにa
0
, a
1
,…,a
n-1
を決定する．
連立1次方程式
の解である．
最小2乗近似：プログラム：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++) { a[i]=0.0; for( j = 0; j < M; j++) { a[i] += y[j]*pow(x[j],(double)i) ; } } /* 係数行列の作成 */ for( i = 0; i <= N; i++) { for( j = 0; j <= i; j++ ) { p[i][j]=0.0; for( k =0; k < M; k++) {
p[i][j] += pow( x[k], (double)(i+j) ); } p[j][i] = p[i][j]; } } /* 連立1次方程式を解く. 結果はaに上書き */ gauss2( p, a, N+1 ); /* 結果の出力 */ 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"); }
最小2乗近似：実行結果 p.123
input_func.dat
0.0 0.2 0.4 0.6 0.8 1.0
2.0 2.12 1.62 2.57 1.53 2.0
output_func.dat
ベクトルxは次の通りです
0.00 0.20 0.40 0.60 0.80 1.00
ベクトルyは次の通りです
2.00 2.12 1.62 2.57 1.53 2.00
最小2乗近似式は
y = 0.38 x^3 - 0.76 x^2 + 0.28 x^1 + 2.00 x^0
0 0.5 1 1.5 2 2.5 3 0 0.2 0.4 0.6 0.8 1 1.2 系列1 系列2ラグランジュ補間 p.124
•点(x
0
,f
0
), (x
1
,f
1
), … , (x
n
,f
n
) が与えられたとき，これらすべての
点(x
i
,f
i
) を通る曲線 y=f(x) を求めて x
0
< x < x
n
の与えられた点
以外の関数値を求めることを補間，内挿 (interpolation)するとういう．
f
k
= f(x
k
), k=0,1,…,n が与えられたとき，等式 P(x
k
) = f
k
, k=0,1,…,n
を満たす多項式 P(x) を f(x) の補間多項式という．
定理6.5 補間条件を満たすn次多項式 P
n
(x) はただ１つに定まる．
Vの行列式：バンデルモンド(Vandermonde)の行列， 解が1つに定まる．
ラグランジュ補間 p.125
•ｎ次のラグランジュ補間多項式，ラグランジュ補間(Lagrange interpolation)
関数近似 確認問題 解答
3.ラグランジュ補間を用いて(1,-1), (2,-1), (3, 7)
を通る2 次曲線を求めよ.
関数近似 確認問題 解答
3.ラグランジュ補間を用いて，を用いて(1,-1), (2,1), (3, 7)
を通る2 次曲線を求めよ.
𝑙𝑙
0
＝
(𝑑𝑑−𝑑𝑑
1
)(𝑑𝑑−𝑑𝑑
2
)
(𝑑𝑑
0
−𝑑𝑑
1
)(𝑑𝑑
0
−𝑑𝑑
2
)
, 𝑙𝑙
1
＝
(𝑑𝑑−𝑑𝑑
0
)(𝑑𝑑−𝑑𝑑
2
)
(𝑑𝑑
1
−𝑑𝑑
0
)(𝑑𝑑
1
−𝑑𝑑
2
)
, 𝑙𝑙
2
＝
(𝑑𝑑−𝑑𝑑
0
)(𝑑𝑑−𝑑𝑑
1
)
(𝑑𝑑
2
−𝑑𝑑
0
)(𝑑𝑑
2
−𝑑𝑑
1
)
𝑙𝑙
0
＝
(𝑑𝑑−2)(𝑑𝑑−3)
2
, 𝑙𝑙
1
＝
−(𝑥𝑥 − 1)(𝑥𝑥 − 3), 𝑙𝑙
2
＝
(𝑑𝑑−1)(𝑑𝑑−2)
2
𝑦𝑦 = 2𝑥𝑥
2
− 4𝑥𝑥 + 1
ラグランジュ補間：プログラム 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;
}
/* ラグランジュ補間 */
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]); } }
