数値解析
2023年度前期 第12週 [7月6日]

静岡大学
工学研究科機械工学専攻
ロボット・計測情報分野
創造科学技術大学院
情報科学専攻三浦 憲二郎講義アウトライン [7月6日]•関数近似と補間•最小2乗近似による関数近似 (復習)
•ラグランジュ補間•形状処理工学の基礎•点列からの曲線の生成

T.Kanai, U.Tokyo

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

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

定理6.1 p.117 （復習）•関数系 φi が１次独立であれば，(6.1)のノルム||・||22を最小にする係数は連立1次方程式

𝑎𝑎0(𝜑𝜑0 , 𝜑𝜑0 )⋯(𝜑𝜑0 , 𝜑𝜑𝑛𝑛−1 )(𝑓𝑓, 𝜑𝜑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.500

0.2

0.4

0.6

0.8

11.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) での値は

ラグランジュ補間：プログラム 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

1086420-2-4-6-8-10系列100.5

11.522.5系列2関数近似 確認問題
1.最小2 乗法を用いて，3 点(0,-1), (1, 1), (2, 7) を近似する直線を求めよ．
2.ラグランジュ補間を用いて(0,-1), (1,1), (2, 7)を通る2 次式を求めよ.
3. 2 次関数y = 2x2 - 1 を0<=x <= 2 の範囲で最良に近似する1 次関数y=a0 + a1 xを求めよ.

関数近似 確認問題 解答
1.最小2 乗法を用いて，3 点(1,-1), (2, 1), (3, 7) を近似する直線を求めよ．
5𝑦𝑦 = 4𝑥𝑥 −3
2.ラグランジュ補間を用いて(0,-1), (1,1), (2, 7)を通る2 次式を求めよ.
𝑦𝑦 = 2𝑥𝑥 2 − 1

3. 2 次関数y = 2x2 - 1 を0<=x <= 2 の範囲で最良に近似する1 次関数y=a0 + a1 xを求めよ.
7
𝑦𝑦 = 4𝑥𝑥 −3

形状処理工学の基礎•形状処理工学 （Computer Aided Geometric Design）
かたちをコンピュータで処理する工学
・表現（データ構造），表示，変形・リバースエンジニアリング
点群から曲面を再構築する．
•ここでは，点群の近似する曲線を生成することを考える．
y=f(x) の形式では，円のように，同じxの値に複数のyの値が対応する関数（多値関数）を表現できない．
そこで, パラメータ t を導入して曲線 C(t)を以下のように定義する．
C(t)=( x(t), y(t) )定理6.4の拡張•m組の与えられた点(x0,y0), (x1,y1), … ,(xm-1,ym-1)とそれらのデータに対応するパラメータ値t0,t1,…,tm-1に対して，それらの点列を内挿する関数C（ｔ）=（x(t),y(t) ）を(gx(t),gy(t))=（∑ai φi(ｔ)，∑bi φi(t))によって近似することを考える．

が最小になるようにa0, a1,…,an-1およびb0, b1,…,bn-1決定する．
2組の連立1次方程式の解である．

点列の補間：プログラム#include <stdio.h>#include <stdlib.h>#include <math.h>#define M 6 /* データのペア数 */#define N 3 /* N次式で近似 *//* ベクトルの入力 */void input_vector2( double *b, char c, FILE *fin, FILE *fout);/* 行列の領域確保 *//* 部分ピボット選択付きガウス消去法 */void gauss2( double a[N+1][N+1], double b[N+1], int n );/* 最小2乗近似呼び出し */void least_square_coef( double x[M], double y[M], double a[N+1], FILE *fout );/* 最小2乗近似 */void least_square( double x[M], double y[M], FILE *fout );/* 曲線の生成 */void curve( double x[M], double y[M], double t[M], FILE *fout);/* 曲線上の点列データ出力 */void curve_shape( double a0[N+1], double b0[N+1], FILE *fp1);
double a[N+1], b[N+1]; /* グローバル変数 */int main(void){
FILE *fin, *fout, *fout1;
double x[M], y[M], t[M];/* ファイルのオープン */if ( (fin = fopen( "curve_fitting.dat", "r")) == NULL ){printf("ファイルが見つかりません : curve_fitting.dat ¥n"); exit(1);}if( (fout = fopen( "output_curve.dat", "w")) == NULL ){printf("ファイルが作成できません : output_curve.dat ¥n"); exit(1);}if( (fout1 = fopen( "output_shape.dat", "w")) == NULL ){printf("ファイルが作成できません : output_shape.dat ¥n");exit(1);}input_vector2( x, 'x', fin, fout ); /* ベクトルxの入出力 */input_vector2( y, 'y', fin, fout ); /* ベクトルyの入出力 */input_vector2( t, 't', fin, fout ); /* パラメータtの入出力 */curve( x, y, t, fout); /* 最小2乗近似 */curve_shape( a, b, fout1);
fclose(fin); fclose(fout); /* ファイルのクローズ */return 0;}

void curve( double x[M], double y[M], double t[M], FILE *fout){
least_square_coef( t, x, a, fout);
least_square_coef( t, y, b, fout);}
void least_square_coef( double x[M], double y[M], double a[N+1],FILE *fout ){double 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 );点列の補間：プログラム その2/* 結果の出力 */fprintf( fout, "最小2乗近似式はy=");for( i = N ; i >= 0 ; i--){if(i==N){
fprintf(fout, "%5.2f x^%d ", a[i],i);}else{
if(a[i]>0){
fprintf(fout, "+ %5.2f x^%d ", a[i],i);}else{
fprintf(fout, "- %5.2f x^%d ",fabs(a[i]),i);}}}
fprintf(fout, "¥n");}/* 部分ピボット選択付きガウス消去法 */void gauss2( double a[N+1][N+1], double b[N+1], int n ){int i, j, k, ip;
double alpha, tmp;
double amax, eps=pow(2.0, -50.0); /* eps = 2^{-50}とする */for( k = 0; k < n-1; k++){/* ピボットの選択 */amax = fabs(a[k][k]); ip = k;for( i = k+1; i < n; i++){if ( fabs(a[i][k]) > amax ){amax = fabs(a[i][k]); ip = i;}}/* 正則性の判定 */if ( amax < eps ) printf("入力した行列は正則ではない!!¥n");/* 行交換 */if ( ip != k){for( j = k; j < n; j++){tmp = a[k][j]; a[k][j]=a[ip][j]; a[ip][j]=tmp;}
tmp = b[k] ; b[k]=b[ip]; b[ip]=tmp;}/* 前進消去 */for( i = k+1; i < n; i++){alpha = - a[i][k]/a[k][k];for( j = k+1; j < n; j++){a[i][j] = a[i][j] + alpha * a[k][j];}b[i] = b[i] + alpha * b[k];}}/* 後退代入 */b[n-1] = b[n-1]/a[n-1][n-1];for( k = n-2; k >= 0; k--){tmp = b[k];for( j = k+1; j <= n; j++){tmp = tmp - a[k][j] * b[j];}b[k] = tmp/a[k][k];}}void curve_shape( double a0[N+1], double b0[N+1], FILE *fp1){
int i,j;double dt0=0.05;
double x0,y0,t0=0.,t1;
for(i=0;i<=20;++i){x0=a0[0];y0=b0[0];t0+=dt0;
t1=t0;
for(j=1;j<=N;++j){x0+=a0[j]*t1;
y0+=b0[j]*t1;
t1*=t0;}fprintf(fp1,"%lf¥t%lf¥n",x0,y0);}}

ラグランジュ補間：実行結果
curve_fitting.dat
0.0 0.5 0.6 0.4 0.8 1.0
2.0 2.57 2.12 1.53 1.62 2.0
0.0 0.2 0.4 0.6 0.8 1.0

32.5
2
1.5

系列1

10.5
000.2

0.4

0.6

0.8

11.2

32.5
2
1.5

系列1

10.5
000.2

0.4

0.60.8

11.2

1.4

連絡先
•電子メイルアドレス miura.kenjiro@shizuoka.ac.jp•ホームページ
https://mc2-lab.com/•電話･ファックス

053-478-1074•授業用ホームページ

https://mc2-lab.com/lecture.html#A

まとめ•関数近似と補間•最小2乗近似による関数近似•ラグランジュ補間•形状処理工学の基礎•点列からの曲線の生成•受講、お疲れ様でした。

T.Kanai, U.Tokyo