6. 関数近似 : 補間と補外 6. ラグランジュ補間法 互いに異なる点 x,x,,x とそれらの点における関数値 f(x ),f(x ),,f(x ) が与えられているとする これらの 点を補間するたかだか - 次の補間多項式 F (x) は, ラグランジュ基底関数 L k (-) (x) を用いて ( ) () F( x) = f( xk) Lk ( x) k= L ( x ) = と書ける これは, 次の多項式 を導入すると ( x x )( x x ) L( x x )( x x ) L( x x ) ( ) k k+ k l ( xk x)( xk x) L( xk xk )( xk xk+ ) L( xk x π ( x) = ( x x )( x x ) L ( x x ) ( ) = ( )( ) L( )( ) L ( ) π x k xk x xk x xk xk xk xk+ xk x であるので,() 式は π (x) がくくり出された形でつぎのようにも表すことができる () F ( x) = π ( x) k= f( xk ) ( x x ) π ( x ) k k ) アルゴリズム 6.. () 式に基づくもの : (x,f ),(x,f ),,(x,f ) := 与えられた 点のデータ fx=0 for k= to p=q= for = to f k the p:=p*(x-x ) q:=q*(x k -x ) ed fx:=fx+f k *p/q ed ed プログラム 6.. 3 4 5 6 7 8 9 /* Lagraga terpolato */ /* gve data: values of (x, f) */ /* terpolato pot xx */ /* terpolated value f(xx) s computed */ #clude <stdo.h> #defe NMAX 0
0 3 4 5 6 7 8 9 0 3 4 5 6 7 8 9 30 3 3 33 34 35 36 37 38 39 40 t ; float x[nmax], f[nmax]; vod ma(vod) t,k; float xx,fx,p,q; prtf("lagraga terpolato "); prtf("umber of data? "); scaf("%d",&); for(=0;<;++) prtf("x[%d], f[%d]? ",+,+); scaf("%g%g", &x[], &f[]); prtf("terpolato pot xx? "); scaf("%g", &xx); fx=0; for(k=0;k<;k++) p=; q=; for(=0;<;++) f(!=k) p *= xx-x[]; q *= x[k]-x[]; fx += f[k]*p/q; prtf(" terpolated value at x=%g: ",xx); prtf("f(%g)=%g ",xx,fx); /* read put data */ /* terpolato */ アルゴリズム 6.. () 式に基づくもの : (x,f ),(x,f ),,(x,f ) := 与えられた 点のデータ for k= to q= for = to f k the q:=q(x k -x ) d k :=q ed ed p=, sum=0 for k= to p :=p(x-x k ) σ:=σ+f k /((x-x k )d k ) ed fx:=pσ
プログラム 6.. 3 4 5 6 7 8 9 0 3 4 5 6 7 8 9 0 3 4 5 6 7 8 9 30 3 3 33 34 35 36 37 38 39 40 4 4 43 /* Lagraga terpolato */ /* gve data: values of (x, f) */ /* terpolato pot xx */ /* terpolated value f(xx) s computed */ #clude <stdo.h> #defe NMAX 0 t ; float x[nmax], f[nmax]; vod ma(vod) t,k; float xx,fx,p,q,sum,pd[nmax]; prtf("lagraga terpolato "); prtf("umber of data? "); scaf("%d",&); for(=0;<;++) prtf("x[%d], f[%d]? ",+,+); scaf("%g%g", &x[], &f[]); prtf("terpolato pot xx? "); scaf("%g", &xx); for(k=0;k<;k++) q=; for(=0;<;++) f(!=k) q *=x[k]-x[]; pd[k]=q; p=; sum=0; for(k=0;k<;k++) p *= xx-x[k]; sum += f[k]/((xx-x[k])*pd[k]); fx=p*sum; prtf(" terpolated value at x=%g: ",xx); prtf("f(%g)=%g ",xx,fx); /* d[k] */ /* terpolato */ 問題 水の動粘性係数 ν は, 常温の範囲では以下のとおりである 3
T [ ] ν[m /s] 0.0.30700E-06 5.0.3900E-06 0.0.00400E-06 5.0 8.9800E-07 30.0 8.00800E-07 ラグランジュ補間により, 水温 8. のときの動粘性係数を求めよ 実行例 プログラム 6.., プログラム 6.. とも実行結果は以下のようになる ------------------------- 実行開始 ----------------------- Lagraga terpolato umber of data? 5 x[], f[]? 0.0.307e-6 x[], f[]? 5.0.39e-6 x[3], f[3]? 0.0.004e-6 x[4], f[4]? 5.0 8.98e-7 x[5], f[5]? 30.0 8.008e-7 terpolato pot xx? 8. terpolated value at x=8.: f(8.)=.04948e-06 ------------------------- おしまい ----------------------- 6. スプライン補間法 互いに異なる点 x 0 < x < < x とそれらの点における関数値 y 0 =f(x 0 ),y =f(x ),,y =f(x ) が与えられているとする 小区間 [x -,x ] ごとに 3 次エルミート補間関数 ( 階導関数の連続性を保証 )s (x) を用いると, x x x x x x x x = + + + + + h h h h s ( x) ( 3y hy ) ( y hy ) ( 3y hy ) ( y hy ) h = x x となる 更に隣接する小区間上の補間関数と 階導関数まで連続に接続する条件 s (x )=s + (x ) を課すると,y に関する方程式 3 3 3 3 y + + y + y + = y + y + y +, =,, L, h h h+ h+ h h h+ h+ を得る 端点 x=x 0,x においては, 階微係数が零となる条件を課すると, 3 3 y + y = y + y h h h h 0 0 3 3 y + y = y + y h h h h となる これらは,y に関する三項方程式となっているので, アルゴリズム 4.3. により y は求まる 点 x がどの区間にあるかを判定すれば,x におけるスプライン補間関数 s (x) の値は定まる プログラム 6.. /* Sple Iterpolato */ /* gve data: umber of tervals: */ 4
3 4 5 6 7 8 9 0 3 4 5 6 7 8 9 0 3 4 5 6 7 8 9 30 3 3 33 34 35 36 37 38 39 40 4 4 43 44 45 46 47 48 49 50 5 5 53 /* x_, y_, =0,,..., */ /* terpolated value f(xx) s computed */ /* for x_0 <= xx <= x_ */ #clude <stdo.h> #clude <math.h> #defe NMAX 50 float yy(float xx,float x0,float x,float y0,float y,float dy0,float dy) float z0,z; z0=(xx-x0)/(x-x0); z=(xx-x)/(x0-x); retur(z*z*((3*y0-(x0-x)*dy0)-z*(*y0-(x0-x)*dy0)) +z0*z0*((3*y-(x-x0)*dy)-z0*(*y-(x-x0)*dy))); vod ma(vod) t,,j,m; float x[nmax],y[nmax],h[nmax],a[nmax],b[nmax],c[nmax],d[nmax]; float x0,dx,xx,yyj; /* read put data */ prtf("sple terpolato "); prtf("umber of tervals? "); scaf("%d",&); for(=0;<=;++) prtf("x[%d], y[%d]? ",,); scaf("%g%g", &x[], &y[]); /* set coeffcets of three-term equatos */ for(=;<=;++) h[]=x[]-x[-]; for(=;<;++) a[]=/h[]+/h[+]; b[]=/h[+]; c[]=/h[]; d[]=3/(h[]*h[])*(y[]-y[-]) +3/(h[+]*h[+])*(y[+]-y[]); /* ed-pot codto at x[0] */ a[0]=/h[]; b[0]=/h[]; c[0]=0; d[0]=3/(h[]*h[])*(y[]-y[0]); /* y"=0 */ /* ed-pot codto at x[] */ a[]=/h[]; b[]=0; c[]=/h[]; d[]=3/(h[]*h[])*(y[]-y[-]); /* y"=0 */ /* a[]=; b[]=0; c[]=0; d[]=(y[]-y[-])/h[]; */ /* solve three-term equatos */ /* forward */ b[0] /= a[0]; d[0] /= a[0]; for(=;<=;++) 5
54 55 56 57 58 59 60 6 6 63 64 65 66 67 68 69 70 7 7 73 74 75 76 77 78 79 a[] -= c[]*b[-]; d[] -= c[]*d[-]; d[] /= a[]; b[] /= a[]; for(=-;>=0;--) d[] -=b[]*d[+]; /* backward */ /* terpolated values */ prtf("tal value ad terval of x for terpolated results, x0, dx? "); scaf("%g%g", &x0,&dx); prtf(" tx t terpolato "); m=(t)floor((x[]-x0)/dx + 0.5); =; /* tal set for terval couter */ for(j=0;j<=m;j++) xx=x0+j*dx; label0:f(xx<x[] ==) yyj= yy(xx,x[-],x[],y[-],y[],d[-],d[]); prtf("%5.6e %5.6e ", xx, yyj); f(xx>=(x[]+0.*dx) && ==) break; else ++; goto label0; 問題 水の動粘性係数 ν の常温での温度依存性は 6. 節 ラグランジュ補間法 の問題で記したとおりである スプライン補間により, 水温が 0~30[ ] のときの動粘性係数を求め, グラフ表示せよ 実行例 ------------------------- 実行開始 ----------------------- Sple terpolato umber of tervals? 4 は小区間数であることに注意 x[0], y[0]? 0.0.307e-6 ( データ点数は+) x[], y[]? 5.0.39e-6 x[], y[]? 0.0.004e-6 x[3], y[3]? 5.0 8.98e-7 x[4], y[4]? 30.0 8.008e-7 tal value ad terval of x for terpolated results, x0, dx? 0.0 0. 初期値 0, 増分 0.のxの値に対する補間値を求める x terpolato.000000e+0.307000e-06.00000e+0.9998e-06.040000e+0.9965e-06.060000e+0.85955e-06.080000e+0.78954e-06.00000e+0.7963e-06.0000e+0.64988e-06 6
.40000e+0.5809e-06.60000e+0.509e-06.80000e+0.4476e-06.00000e+0.3786e-06 : : : :.900000e+0 8.8459e-07.90000e+0 8.4894e-07.940000e+0 8.363e-07.960000e+0 8.078395e-07.980000e+0 8.04390e-07 3.000000e+0 8.008000e-07 -------------------------おしまい----------------------- 動粘性係数の温度依存性はつぎのグラフに示される 動粘性係数.40E-06.30E-06 動粘性係数 [m^/sec].0e-06.0e-06.00e-06 9.00E-07 8.00E-07 7.00E-07 0 5 0 5 30 温度 [ ] 参考文献 葉子 : 数値計算の基礎 解法と誤差, コロナ社 (007) 森口繁一, 伊理正夫, 武市正人編 :C による算法痛論, 東京大学出版会 (000) Heath, Mchael T.: Scetfc Computg, A Itroductory Survey, McGraw-Hll(00) 7