290 likes | 428 Views
数値解析シラバス. C言語環境と数値解析の概要(1回) 方程式の根(1回) 連立一次方程式(2回) 補間と近似(2回) 数値積分(1回) 中間試験(1回) 常微分方程式(1回) 偏微分方程式(2回) 数値計算の応用(3回) 期末試験(1回). 最小二乗法. 実験データなど誤差を含む点を近似する. 補間 データ点上を必ず通過する近似. 近似(最小2乗法) データ点に近い点を通る近似. y. y. x. x. 関数の近似法. 最小二乗法. y. y = a 0 + a 1 x + a 2 x 2. D y 3. D y 0.
E N D
数値解析シラバス • C言語環境と数値解析の概要(1回) • 方程式の根(1回) • 連立一次方程式(2回) • 補間と近似(2回) • 数値積分(1回) • 中間試験(1回) • 常微分方程式(1回) • 偏微分方程式(2回) • 数値計算の応用(3回) • 期末試験(1回)
最小二乗法 実験データなど誤差を含む点を近似する
補間 データ点上を必ず通過する近似 近似(最小2乗法) データ点に近い点を通る近似 y y x x 関数の近似法
最小二乗法 y y = a0 + a1x + a2x2 Dy3 Dy0 Dy1 Dy2 x 最小になるように a0, a1,a2を決定 Dy02 + Dy12 + Dy22 + Dy32
近似(最小2乗法) 1次関数で 一次関数の例 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 0 この和を最小にする a0, a1
2 i=0 SDyi2= (a0+a1x0–y0)2+(a0+a1x1–y1)2 +(a0+a1x2–y2)2 a0,a1で偏微分 m i=0 ∂ ∂a0 S Dyi2= 2(a0+a1x0–y0) +2(a0+a1x1–y1)+2(a0+a1x2–y2) = 0 m i=0 ∂ ∂a1 S Dyi2= 2(a0+a1x0–y0)x0 +2(a0+a1x1–y1)x1+2(a0+a1x2–y2)x2 = 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)x1a1 (x0y0+x1y1+x2y2)
二次関数の例 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)
近似(最小2乗法) 1次関数で 紙とペン y = a + bx y 1 x 0 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とおく。 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 (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 y 1 3a + 3b – 2 = 0 b = 1/2 傾き:1/2 交点:1/6 a= 1/6 3a+ 5b -3=0 x 0 1 2
近似(最小2乗法) 1次関数で 紙とペン y = a + bx y 1 x 0 1 2
紙とペン aで偏微分→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 (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 y 傾き:1/2 交点: – 1/6 1 3a + 3b – 1 = 0 b = 1/2 3a+ 5b – 2=0 a= – 1/6 x 0 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; プログラミング
m k = 0 S xki+j 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)); for(i=0;i<=M;i++) for(k=0;k<N;k++) a[i][M+1] += y[k]*pow(x[k],(double)i); if(jordan()==1) return 1; for(i=0;i<=M;i++) printf("A%2d = %7.3lf\n", i, a[i][M+1]); return 0; } m k = 0 S xkiyk 1のときは,エラーが起きている
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]; } 連立方程式のプログラミング エラーが起こると 1を返す.
テキストのままで実行 テキストと異なる A0 = 0.002 A1 = 0.772 A2 = 0.395 A3 =‐0.079 A4 = 0.002 しかも。エラーがでる。
エクセルで確認 答えは正しい!
エラーの原因 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)); for(i=0;i<=M;i++) for(k=0;k<N;k++) a[i][M+1] += y[k]*pow(x[k],(double)i); pow(x,0)で エラー発生
エラーの回避法 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)); } for(i=0;i<=M;i++) for(k=0;k<N;k++) if((i)<EPS){ a[i][M+1] += y[k]; }else{ a[i][M+1] += y[k]*pow(x[k],(double)i); } pow(x,0)で powを1とする. pow(x,0)で powを1とする.
実行例 エラーがなくなった!
演習問題3.2 最小二乗法で直線近似 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}; 変更 変更
演習問題3.4 最小二乗法で3次方程式近似 Cプログラムの後 エクセルで確認