数値解析
2022年度前期 第11週 [6月29日]

静岡大学
創造科学技術大学院情報科学専攻
工学部機械工学科ロボット・計測情報分野
三浦 憲二郎講義アウトライン [6月29日]•関数近似と補間•最小2乗近似による関数近似•ラグランジュ補間T.Kanai, U.Tokyo

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

関数近似•閉区間[a,b]で定義された関数f(x)をｇ（ｘ）=∑aiφi(x)で近似する．関数系 φi(x) は[a, b]上で連続かつ１次独立関数が1次独立とは，関数の差のノルムの二乗が最小になるように係数 ai を定める．

関数近似 確認問題
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=-52.次の関数の最小値を求めよ．
𝑑𝑑𝑑𝑑
𝑑𝑑𝑑𝑑
= 2𝑥𝑥 + 𝑦𝑦 + 1 = 0,= 𝑥𝑥 + 2𝑦𝑦 = 0
𝑑𝑑𝑑𝑑

21𝑥𝑥 = , 𝑦𝑦 = − ,33𝑑𝑑𝑦𝑦3.最小2 乗法を用いて，3 点(1,-1), (2, 1), (3, 7) を近似する
直線を求めよ．
17𝑦𝑦 = 4𝑎𝑎 −3

定理6.1 p.117•関数系 φi が１次独立であれば，(6.1)のノルム||・||22を最小にする係数は連立1次方程式(𝜑𝜑0 , 𝜑𝜑0 )⋮(𝜑𝜑𝑛𝑛−1 , 𝜑𝜑0 )𝑎𝑎0⋯(𝜑𝜑0 , 𝜑𝜑𝑛𝑛−1 )(𝑓𝑓, 𝜑𝜑0 )⋮ =⋱⋮⋮⋯ (𝜑𝜑𝑛𝑛−1 , 𝜑𝜑𝑛𝑛−1 ) 𝑎𝑎𝑛𝑛−1 (𝑓𝑓, 𝜑𝜑𝑛𝑛−1 )の解として一意に決定される．ここで，[a, b] 上の連続関数f(x)とg(x)に対して

定理6.1 p.117•(証明) (6.1)よりが最小になるには極値の必要条件を満たさなければならない．

連立１次方程式(6.2）を解くことと同値．

グラム行列式 p.117•要素 x0, x1, …,xn-1 が線形独立であることと，以下のグラム行列式が0でないことと同値．（6.2）はただ１つの解を持つ．

定理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乗近似：プログラム：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

32.5

2

output_func.dat
ベクトルxは次の通りです
0.00 0.20 0.40 0.60 0.80 1.00
1.5ベクトルyは次の通りです
2.00 2.12 1.62 2.57 1.53 2.001
最小2乗近似式は
y = 0.38 x^3 - 0.76 x^2 + 0.28 x^1 + 2.00 x^0

系列1系列20.5

000.20.4

0.60.811.2

ラグランジュ補間 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)の行列， 解が1つに定まる．

ラグランジュ補間 p.125•ｎ次のラグランジュ補間多項式，ラグランジュ補間(Lagrange interpolation)

基本多項式 li(x) の 点xk (0≦k≦n) での値は関数近似 確認問題 解答
3.ラグランジュ補間を用いて(1,-1), (2,-1), (3, 7)を通る2 次曲線を求めよ.

関数近似 確認問題 解答
3.ラグランジュ補間を用いて，を用いて(1,-1), (2,1), (3, 7)を通る2 次曲線を求めよ.
(𝑥𝑥−𝑥𝑥1 )(𝑥𝑥−𝑥𝑥2 )(𝑥𝑥−𝑥𝑥0 )(𝑥𝑥−𝑥𝑥2 )(𝑥𝑥−𝑥𝑥0 )(𝑥𝑥−𝑥𝑥1 )𝑙𝑙0 ＝, 𝑙𝑙 ＝, 𝑙𝑙 ＝(𝑥𝑥0 −𝑥𝑥1 )(𝑥𝑥0 −𝑥𝑥2 ) 1(𝑥𝑥1 −𝑥𝑥0 )(𝑥𝑥1 −𝑥𝑥2 ) 2(𝑥𝑥2 −𝑥𝑥0 )(𝑥𝑥2 −𝑥𝑥1 )(𝑥𝑥−2)(𝑥𝑥−3)(𝑥𝑥−1)(𝑥𝑥−2), 𝑙𝑙1 ＝−(𝑥𝑥 − 1)(𝑥𝑥 − 3), 𝑙𝑙2 ＝
22

𝑙𝑙0 ＝

𝑦𝑦 = 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 );/* ラグランジュ補間 */double lagrange( double x[N+1], double y[N+1], double xi ){int i, k;double pn = 0.0, li;

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;}/* 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]);}}ラグランジュ補間：実行結果
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.1 1.6 2.6 1.5 2.7 0.67 3.5 0.94 2.0

10864
20-2-4-6-8-10

系列100.5

11.522.5系列2関数近似 確認問題
４. 2 次関数y = 2x2 - 4x + 1 を1 <=x <= 3 の範囲で最良に近似する1 次関数y=a0 + a1 xを求めよ.

関数近似 確認問題 解答
４. 2 次関数y = 2x2 - 4x + 1 を1 <=x <= 3 の範囲で最良に近似する1 次関数y=a0 + a1 xを求めよ.
24426310𝑎𝑎0
3=𝑎𝑎1 28319𝑦𝑦 = 4𝑥𝑥 −3まとめ•関数近似と補間•最小2乗近似による関数近似•ラグランジュ補間T.Kanai, U.Tokyo