数値解析シラバス Ｃ言語環境と数値解析の概要（１回） 方程式の根（１回） 連立一次方程式（２回） 補間と近似（２回） 数値積分（１回） 中間試験（１回） 常微分方程式（１回） 偏微分方程式（２回） 数値計算の応用（３回） 期末試験（１回）

最小二乗法 実験データなど誤差を含む点を近似する

関数の近似法 補間 　データ点上を必ず通過する近似 近似（最小２乗法） 　データ点に近い点を通る近似 x y x y

最小二乗法 y = a0 + a1x + a2x2 Dy3 Dy0 Dy1 Dy2 Dy02 + Dy12 + Dy22 + Dy32 最小になるように a0, a1,a2を決定 Dy02 + Dy12 + Dy22 + Dy32

一次関数の例 y = a0 + a1 x Dy02 = (a0 + a1x0 - y0)2 Dy12 = (a0 + a1x1 - y1)2 近似（最小２乗法） 　１次関数で y = a0 + a1 x y Dy02 = (a0 + a1x0 - y0)2 (x1, y1) Dy12 = (a0 + a1x1 - y1)2 (x0, y0) (x2, y2) Dy22 = (a0 + a1x2 - y2)2 y = a0 + a1 x x この和を最小にする a0, a1

SDyi2 = (a0+a1x0–y0)2+(a0+a1x1–y1)2 +(a0+a1x2–y2)2 m i=0 ∂ ∂a0 = 0 S Dyi2 = 2(a0+a1x0–y0)x0 +2(a0+a1x1–y1)x1+2(a0+a1x2–y2)x2 m i=0 ∂ ∂a1 = 0

= (a0+a1x0–y0)+(a0+a1x1–y1)+(a0+a1x2–y2) =(3)a0+(x0+x1+x2)a1–(y0+y1+y2)=0 (a0+a1x0–y0)x0+(a0+a1x1–y1)x1+(a0+a1x2–y2)x0 =(3)a0x0+(x0+x1+x2)a1x1–(x0y0+x1y1+x2y2)=0 (3)a0+(x0+x1+x2)a1=(y0+y1+y2) (3)a0x0+(x0+x1+x2)a1x1=(x0y0+x1y1+x2y2) (3) (x0+x1+x2) a0 (y0+y1+y2) = (3)x0 (x0+x1+x2)x1 a1 (x0y0+x1y1+x2y2)

二次関数の例 ∂ S Dyi2 = S 2xij (a0 + a1xi + a2xi2 - yi) = 0 ∂aj m i=0 m i=0 ∂ ∂aj S Dyi2 = S 2xij (a0 + a1xi + a2xi2 - yi) = 0 追加 (x0j +･･･+ xmj )a0 + (x0j+1 +･･･+ xmj+1 ) a1 + (x0j+2 +･･･+ xmj+2 ) a2 = (x0j y0 +･･･+ xmj ym ) j = 0, 1, 2

(x00 +･･･+ xm0)a0+(x01 +･･･+ xm1)a1+(x02 +･･･+ xm2) a2 = (x00 y0 +･･･+ xm0 ym ) (x01+･･･+ xm1)a0+(x02 +･･･+ xm2)a1+(x03 +･･･+ xm3) a2 = (x01 y0 +･･･+ xm1 ym ) (x02+･･･+ xm2)a0+(x03 +･･･+ xm3)a1+(x04 +･･･+ xm4) a2 = (x02 y0 +･･･+ xm2 ym )

j = 0, 1, 2 (x00 +･･･+ xm0) (x01 +･･･+ xm1) (x02 +･･･+ xm2) a0 (x01+･･･+ xm1) (x02 +･･･+ xm2) (x03 +･･･+ xm3) a1 (x02+･･･+ xm2) (x03 +･･･+ xm3) (x04 +･･･+ xm4) a2 i = 0, 1, 2 (x00 y0 +･･･+ xm0 ym ) m k = 0 m k = 0 = (x01 y0 +･･･+ xm1 ym ) S xkiyk S xki+j (x02 y0 +･･･+ xm2 ym )

紙とペン 近似（最小２乗法） 　１次関数で y = a + bx x y 1 2 x y 1 2

紙とペン x y 1 2 Dy02 = (a + bx0 - y0)2 Dy12 = (a + bx1 - y1)2 1 2 Dy02 = (a + bx0 - y0)2 Dy12 = (a + bx1 - y1)2 Dy22 = (a + bx2 - y2)2 Dy02 + Dy12 +Dy22 = (a + bx0 - y0)2 +(a + bx1 - y1)2 +(a + bx2 - y2)2 ⇒a,bで偏微分

紙とペン aで偏微分→0とおく。 bで偏微分→0とおく。 2(a + bx0 - y0)+2(a + bx1 - y1) +2(a + bx2 - y2)=0 bで偏微分→0とおく。 2(a + bx0 - y0) x0+2(a + bx1 - y1) x1 +2(a + bx2 - y2) x2=0 x y 1 2 (a)+(a + b - 1) +(a + 2b - 1)=0 (a) + (a + b - 1)+2(a + 2b – 1)=0

紙とペン (a)+(a + b - 1) +(a + 2b - 1)=0 (a) + (a + b - 1)+2(a + 2b – 1)=0 x y 1 2 (a)+(a + b - 1) +(a + 2b - 1)=0 (a) + (a + b - 1)+2(a + 2b – 1)=0 3a + 3b – 2 = 0 3a + 5b -3=0 x y 1 2 b = 1/2 a = 1/6 傾き：1/2 交点：1/6

紙とペン 近似（最小２乗法） 　１次関数で y = a + bx x y 1 2 y 1 x 1 2

紙とペン aで偏微分→0とおく。 bで偏微分→0とおく。 2(a + bx0 - y0)+2(a + bx1 - y1) +2(a + bx2 - y2)=0 bで偏微分→0とおく。 2(a + bx0 - y0) x0+2(a + bx1 - y1) x1 +2(a + bx2 - y2) x2=0 x y 1 2 (a)+(a + b) +(a + 2b - 1)=0 (a + b )+2(a + 2b – 1)=0

紙とペン (a)+(a + b ) +(a + 2b - 1)=0 (a + b )+2(a + 2b – 1)=0 x y 1 2 (a)+(a + b ) +(a + 2b - 1)=0 (a + b )+2(a + 2b – 1)=0 3a + 3b – 1 = 0 3a + 5b – 2=0 y 傾き：1/2 交点： – 1/6 1 b = 1/2 a = – 1/6 x 1 2

プログラミング #include <stdio.h> #include <math.h> #define N 6 #define M 4 #define EPS .00001 double a[M+1][M+2]; int jordan( void ); int main( ) { double x[N]={0.0, 1.0, 2.0, 3.0, 3.1, 5.0}; double y[N]={0.0, 1.1, 2.5, 4.0, 4.1, 5.0}; int i, j, k; for(i=0;i<=M;i++) for(j=0;j<=M+1;j++) a[i][j]=0.0; プログラミング

S xki+j S xkiyk m for(i=0;i<=M;i++) k = 0 for(j=0;j<=M;j++) for(k=0;k<N;k++) a[i][j] += pow(x[k],(double)(i+j)); a[i][M+1] += y[k]*pow(x[k],(double)i); if(jordan()==1) return 1; printf("A%2d = %7.3lf\n", i, a[i][M+1]); return 0; } S xki+j m k = 0 S xkiyk １のときは，エラーが起きている

連立方程式のプログラミング int jordan( void ){ double pivot, del; int i, j, k; for(i=0;i<=M;i++) { pivot = a[i][i]; if( fabs(pivot) < EPS ){ printf("ピボット数が許容範囲以下\n"); return 1; } for( j=i; j <= M+1; j++) a[i][j] /= pivot; for(k=0; k<=M;k++){ del = a[k][i]; for( j=0;j<=M+1; j++){ if( k!=i) a[k][j] -= del*a[i][j]; 連立方程式のプログラミング エラーが起こると １を返す．

テキストのままで実行 A0 = 0.002 A1 = 0.772 A2 = 0.395 A3 =‐0.079 A4 = 0.002 テキストと異なる A0 = 0.002 A1 =　 0.772 A2 = 0.395 A3 =‐0.079 A4 = 　0.002 しかも。エラーがでる。

エクセルで確認 答えは正しい！

エラーの原因 pow(x,0)で エラー発生 pow(x,0)で エラー発生 for(i=0;i<=M;i++) for(j=0;j<=M;j++) for(k=0;k<N;k++) a[i][j] += pow(x[k],(double)(i+j)); a[i][M+1] += y[k]*pow(x[k],(double)i); pow(x,0)で エラー発生

エラーの回避法 pow(x,0)で powを1とする． pow(x,0)で powを1とする． for(i=0;i<=M;i++) for(j=0;j<=M;j++) for(k=0;k<N;k++) if((i+j)<EPS){ a[i][j] += 1.0; }else{ a[i][j] += pow(x[k],(double)(i+j)); } if((i)<EPS){ a[i][M+1] += y[k]; a[i][M+1] += y[k]*pow(x[k],(double)i); pow(x,0)で powを1とする． pow(x,0)で powを1とする．

実行例 エラーがなくなった！

演習問題３．２ 最小二乗法で直線近似 x 0 1 2 3 4 y 1 2 1 0 4

解答 #include <stdio.h> #include <math.h> #define N 5 #define M 1 #define EPS .00001 double a[M+1][M+2]; int jordan( void ); int main( ) { double x[N]={0.0, 1.0, 2.0, 3.0, 4.0}; double y[N]={1.0, 2.0, 1.0, 0.0, 4.0}; 変更 変更

演習問題３．４ 最小二乗法で３次方程式近似 Ｃプログラムの後 エクセルで確認 x -4 -2 -1 1 3 4 6 y -35.1 15.1 1 3 4 6 y -35.1 15.1 15.9 8.9 0.1 21.1 135.0 Ｃプログラムの後 エクセルで確認

エクセルで結果を確認