Matlab數(shù)值算法實(shí)現(xiàn)

1、Ch10. 數(shù)值算法實(shí)現(xiàn)1. 線性方程組解法 1、三角形線性方程組解法 以上三角形線性方程組 為例。 回代：1 %文件 uptri.mfunction u = uptri ( a, b )n= size(a,1) ;x(n)= b(n) / a(n, n) ;for i = n-1:-1:1 s=0 ; for j = i+1: n s=s+a(i, j) * x( j ) ; end x(i) = ( b(i) s) / a(i, i) ;endu = x ;求解22、順序Gauss消去法 （1）消去過(guò)程：第 k 步，計(jì)算（2）回代過(guò)程： 3%文件 gauss.mfunction u = g

2、auss (a, b)n = length (b) ;for k=1: n 1 for i = k+1 : n mult = a ( i, k) / a (k, k) ; for j =k +1: n % if abs ( a( k, k) ) 1e6 a (i, j) = a (i, j) mult * a(k, j) ; % else % disp (順序Gauss消去法失敗); % pause ； % exit ； % end end可去掉%，判斷主元是否為04 b (i) = b (i) mult * b (k) ; endendx(n)= b(n) / a(n, n) ;for i

3、= n-1:-1:1 s=0 ; for j = i+1: n s=s+a(i, j) * x( j ) ; end x(i) = ( b(i) s) / a(i, i) ;endu = x ;回代5例：% 主文件main.ma=6, -2, 2, 4; 12, -8, 6, 10; 3,-13, 9, 3; -6, 4, 1, -18 ;b=16, 26, -19, -34;x= gauss (a, b);disp ( 方程組解為: );x則有: main方程組解為:x= 3 1 -2 163、Jacobi迭代法 7% jacobi.mfunction y = jacobi ( a, b,

4、x0)D = diag ( diag (a) ) ;U = - triu (a, 1) ;L = - tril (a, -1) ;B = D ( L+U) ;f = D b ;y = B* x0 + f ; n=1;while norm (y x0 ) = 1.0e 6 x0 = y ; y = B * x0 + f ; n = n +1;endyn8例: P249 例7.21 a =10, -1, 0; -1, 10, -2; 0, -2, 10; b =9; 7; 6; jacobi ( a, b, 0; 0; 0 )y = 0.9958 0.9579 0.7916n = 11解向量迭代次

5、數(shù)92. 方程求根1、二分法% erfen.mfunction y = erfen (fun,a, b, esp )if nargin 4 esp =1e 4 ; endif feval (fun, a) * feval (fun, b) esp if feval ( fun, a) * feval (fun, c) 0 b = c ; c = ( a+b) / 2 ; elseif feval ( fun, c) * feval (fun, b) erfen ( fc , 0, 10, 0.05 )n = 56ans = 1.6180先編寫(xiě)函數(shù)122、不動(dòng)點(diǎn)迭代法 % iterate.mfu

6、nction y = iterate (x)x1 = g (x) ;n =1 ;while ( abs ( x1 x ) = 1.0e 6 ) & ( n iterate ( 3 )x1 = 3.7331n = 22初值不同,迭代步數(shù)會(huì)不同迭代函數(shù)143、牛頓迭代法 % newton.mfunction y = newton( x0 )x1 = x0 fc ( x0 ) / df ( x0 ) ; n = 1;while (abs(x1 x0) = 1.0e 6 )& ( n newton (0) newton ( 10 ) x1= x1= -0.4590 3.7331n = n = 6 12

7、比不動(dòng)點(diǎn)迭代快163. 插值方法 以牛頓插值為例。 已知函數(shù) 在互異點(diǎn) 上的值為 。 n次Newton插值多項(xiàng)式 為：17定義 設(shè)函數(shù) 在互異點(diǎn) 上的值為 。定義：1）一階均差為2）二階均差為3）遞推下去，k 階均差為18例： 已知 f (x) = sqrt (2x) + ln (x) 節(jié)點(diǎn) x1 = 0.5 , x2 = 1 , x3 = 1.5 , x4 = 2.0 求 Newton 插值多項(xiàng)式，并計(jì)算其在 x = 0.75 , x = 1.75處的值。比較與原函數(shù)的誤差。19% f.mfunction y = f (x)y = sqrt (2*x) + log (x) ;% j0.mfu

8、nction r = j0 (x)r = f (x) ;% j1.mfunction r = j1 (x1, x2)r = (j0 (x1)-j0(x2) / (x1-x2) ;原函數(shù)0 階均差 (原函數(shù))1 階均差20% j2.mfunction r = j2 (x1, x2, x3)r = (j1 (x1, x2) - j1 (x2, x3) / (x1-x3) ;% j3.mfunction r = j3 (x1, x2, x3, x4)r = (j2 (x1, x2, x3) - j2 (x2, x3, x4) / (x1-x4) ;2 階均差3 階均差21% newton.mfunc

9、tion r = newton ( t )x = 0.5, 1.0, 1.5, 2.0 ;r = j0 (x(1)+ j1 (x(1), x(2).*(t-x(1)+j2 (x(1), x(2), x(3) .*(t-x(1) .*(t-x(2)+ j3 (x(1), x(2), x(3), x(4) .*(t-x(1) .*(t-x(2) .*(t-x(3) ;插值函數(shù)22% main.mdisp ( x=0.75處值與誤差);t =0.75; l =newton(t)e=abs(f(t)-l)disp ( x=1.75處值與誤差);t =1.75; l =newton(t)e=abs(f(

10、t)-l)row =0.5: 0.1: 2.0;y= f(row);N= newton(row);plot (row, y, b-,row, N, r-);legend (原函數(shù), 插值函數(shù));23 mainx=0.75處值與誤差l = 0.9221e = 0.0150 x=1.75處值與誤差l = 2.4228e = 0.0077241、 多項(xiàng)式擬合 已知 ( xi , yi ) ( i = 1, , m) polyfit(x, y, n)x= x1 , x2 , , xm y= y1 , y2 , , ym 多項(xiàng)式次數(shù)25例: 用三次多項(xiàng)式在 0, 5 上擬合 ex x = 0: 0.1:

11、 5; y = exp ( x ); p = polyfit ( x, y, 3 ); s = polyval (p, x ); plot ( x, y, b*, x, s, “r- ) legend ( 原曲線, 擬合曲線 ) axis ( 0, 5, 0, 52 )擬合函數(shù)262、一般情形對(duì)應(yīng)于解超定方程組 用矩陣除法實(shí)現(xiàn)例：用 y = a + bx2 擬合給定的一組數(shù)據(jù)： 27解：使 最小。對(duì)應(yīng)的超定方程組為：即Matlab實(shí)現(xiàn)： u = A y ( 見(jiàn) P228 例7.8 ) 記為 Au = y285. 數(shù)值積分 復(fù)化梯形公式的實(shí)現(xiàn): trapz (y) * h 即y = f(x0),

12、 , f(xn) 29例: P233 例7.10%fun.mfunction y = fun (t)y = exp (-0.5*t).*sin (t+pi/6) ; d= pi /1000 ; t= 0: d: 3*pi; y= fun (t) ; z= trapz (y) * dz = 0.988被積函數(shù)305. 常微分方程數(shù)值解法 以Euler法為例。例:取步長(zhǎng) h = 0.231% f.mfunction r = f (x, y)r = y- x+ 1 ;% orig.mfunction r = orig (x, y)r = x+ exp (x) ;右端函數(shù)準(zhǔn)確解32% euler.m

13、Euler方法的函數(shù)function X, Y = euler (a, b, h, c) % c= y0n = round (b-a)/h) + 1 ; % 計(jì)算點(diǎn)的個(gè)數(shù)x = zeros (n, 1) ; % 設(shè)定 x 維數(shù)y = zeros (n, 1) ; % 設(shè)定 y 維數(shù)x(1) = a ;y(1) = c ; % 初始條件str = sprintf (x0=%g, y0=%g n, x(1), y(1);disp (str); % 顯示初始條件33% 對(duì)每個(gè)點(diǎn)用Euler公式計(jì)算for i=1: n-1 y(i+1)=y(i)+h*f(x(i), y(i) ; x(i+1)=x(i)+h ; str= sprintf (%d x= %g y= %g e= %g, i , x(i+1), y(i+1), abs(y(i+1) - orig(x(i+1) ; disp (str); % 顯示計(jì)算結(jié)果與誤差endX = x ; Y = y ; % 返回計(jì)算結(jié)果34% 主文件 main.mclear ; % 清除內(nèi)存變量a=0; b=5; h=0.2; c