1 / 29

数値解析シラバス

数値解析シラバス. 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.

tacita
Download Presentation

数値解析シラバス

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript


  1. 数値解析シラバス • C言語環境と数値解析の概要(1回) • 方程式の根(1回) • 連立一次方程式(2回) • 補間と近似(2回) • 数値積分(1回) • 中間試験(1回) • 常微分方程式(1回) • 偏微分方程式(2回) • 数値計算の応用(3回) • 期末試験(1回)

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

  3. 補間  データ点上を必ず通過する近似 近似(最小2乗法)  データ点に近い点を通る近似 y y x x 関数の近似法

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

  5. 近似(最小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

  6. 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

  7. (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)

  8. 二次関数の例 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

  9. (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)

  10. 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)

  11. 近似(最小2乗法)  1次関数で 紙とペン y = a + bx y 1 x 0 1 2

  12. 紙とペン 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で偏微分

  13. 紙とペン 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

  14. 紙とペン (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

  15. 近似(最小2乗法)  1次関数で 紙とペン y = a + bx y 1 x 0 1 2

  16. 紙とペン 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

  17. 紙とペン (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

  18. #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; プログラミング

  19. 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のときは,エラーが起きている

  20. 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を返す.

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

  22. エクセルで確認 答えは正しい!

  23. エラーの原因 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)で エラー発生

  24. エラーの回避法 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とする.

  25. 実行例 エラーがなくなった!

  26. 演習問題3.2 最小二乗法で直線近似 x 0 1 2 3 4 y 1 2 1 0 4

  27. 解答 #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}; 変更 変更

  28. 演習問題3.4 最小二乗法で3次方程式近似 Cプログラムの後 エクセルで確認

  29. エクセルで結果を確認

More Related