Matlab数值算法实现.ppt

单击此处编辑母版标题样式,,单击此处编辑母版文本样式,,第二级,,第三级,,第四级,,第五级,,,,*,Ch10. 数值算法实现,§,1. 线性方程组解法,,,1、三角形线性方程组解法,,,以上三角形线性方程组 为例回代：,,,1,,%文件 uptri.m,,function 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) ;,,end,,u = x ;,求解,,,2,,2、,顺序,Gauss,消去法,,（1）,消去过程：,,第 k 步，计算,,,,,,（2）回代过程：,,,,,,,3,,%文件 gauss.m,,function u = gauss (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) ),>,1e–6,,a (i, j) = a (i, j) – mult * a(k, j) ;,,% else,,% disp (‘,顺序,Gauss,消去法失败’);,,,% pause ；,,% exit ；,,% end,,end,可去掉%，,,判断主元是否为0,,4,,b (i) = b (i),–,mult * b (k) ;,,end,,end,,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) ;,,end,,u = x ;,回代,,5,,例：,,% 主文件main.m,,a=[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 1,,6,,3、Jacobi,迭代法,,,,,,,,,,,7,,% jacobi.m,function y = jacobi ( a, b, 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;,,end,,y,,n,,8,,例: 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.7916,,n =,,11,解向量,迭代次数,,9,,§,2. 方程求根,1,、二分法,,% erfen.m,,function y = erfen (fun,a, b, esp ),,if nargin < 4 esp =1e –4 ; end,,if feval (fun, a) * feval (fun, b) < 0,,n = 1 ; c = (a+ b) / 2 ;,,while (b-a) > esp,,if feval ( fun, a) * feval (fun, c) <0,,b = c ; c = ( a+b) / 2 ;,,elseif feval ( fun, c) * feval (fun, b) <0,,a = c ; c = ( a+b) / 2 ;,,else y = c ; esp = 10000 ;,,end,,n= n+1 ;,,end,,10,,,y = c ;,,elseif feval ( fun, a) == 0,,y = a ;,,elseif feval ( fun, b) == 0,,y = b ;,,else,,disp (‘there may not be a root in the interval’);,,end,,n,,11,,例: P258 例7.32,,% fc.m,,function y = fc (x),,Y = x^2 –x –1 ;,,,>> erfen ( ‘fc’ , 0, 10, 0.05 ),,n =,,56,,ans =,,1.6180,,先编写函数,,12,,2、不动点,迭代法,,,% iterate.m,,function y = iterate (x),,x1 = g (x) ;,,n =1 ;,,while ( abs ( x1– x ) >= 1.0e –6 ) & ( n<=1000 ),,x = x1 ;,,x1 = g (x) ; n = n+1 ;,,end,,x1,,n,,13,,例：P259 例7.33,x = ln (3x,2,),迭代函数,(x) = ln (3x,2,),,% g.m,,function y = g(x),,y=log (3*x.^2) ;,,,>> iterate ( 3 ),,x1 =,,3.7331,,n =,,22,初值不同,,,迭代步数会不同,迭代函数,,14,,3、牛顿,迭代法,,% newton.m,,function y = newton( x0 ),,x1 = x0 – fc ( x0 ) / df ( x0 ) ;,,n = 1;,,while (abs(x1 –x0) >= 1.0e –6 )& ( n <= 100000000),,x0 = x1 ;,,x1 = x0 – fc ( x0 ) / df ( x0 ) ; n = n +1;,,end,,x1,,n,,fc 为原始函数,,df 为导函数,,15,,例： P260 例7.34,,% fc.m,,function y = fc (x),,y = 3* x.^2 –exp (x);,,% df.m,,function y = df (x),,y = 6* x –exp (x);,,,>> newton (0) >> newton ( 10 ),,x1= x1=,,-0.4590 3.7331,,n = n =,,6 12,比不动点,迭代快,,16,,§,3. 插值方法,以牛顿插值为例。

已知函数 在互异点 上的值为,,n次,Newton,插值多项式,为：,,17,,定义,,设函数 在互异点 上的值为,,定义：,,1）,一阶均差,为,,,,2）,二阶均差,为,,,,3）递推下去，k,阶均差,为,,,,18,,例：,,已知 f (x) = sqrt (2x) + ln (x),,节点 x,1,= 0.5 , x,2,= 1 , x,3,= 1.5 , x,4,= 2.0,,求 Newton,插值多项式，并计算其在,,x = 0.75 , x = 1.75处的值比较与原函数的误差19,,% f.m,,function y = f (x),,y = sqrt (2*x) + log (x) ;,,,% j0.m,,function r = j0 (x),,r = f (x) ;,,,% j1.m,,function r = j1 (x1, x2),,r = (j0 (x1)-j0(x2)) / (x1-x2) ;,,原函数,0 阶均差 (原函数),1 阶均差,,20,,% j2.m,,function r = j2 (x1, x2, x3),,r = (j1 (x1, x2) - j1 (x2, x3)) / (x1-x3) ;,,,% j3.m,,function r = j3 (x1, x2, x3, x4),,r = (j2 (x1, x2, x3) - j2 (x2, x3, x4)) / (x1-x4) ;,,,2 阶均差,3 阶均差,,21,,% newton.m,,function 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)) ;,,插值函数,,22,,% main.m,,disp ( ‘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(t)-l),,row =0.5: 0.1: 2.0;,,y= f(row);,,N= newton(row);,,plot (row, y, ‘b-’,row, N, ‘r--’);,,legend (‘原函数’, ‘插值函数’);,,23,,>> main,,x=0.75处值与误差,,l =,,0.9221,,e =,,0.0150,,x=1.75处值与误差,,l =,,2.4228,,e =,,0.0077,,,24,,1、 多项式拟合,已知 (,x,i,, y,i,),,( i = 1, …, m),,,polyfit(x, y, n),,,x=[ x,1,, x,2,, … , x,m,],y=[ y,1,, y,2,, … , y,m,],多项式次数,,25,,例: 用三次多项式在 [0, 5] 上拟合,e,x,>> x = 0: 0.1: 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 ] ),,拟合函数,,26,,2、一般情形,对应于解超定方程组,,—用矩阵除法实现,,,例：用 y = a + bx,2,拟合给定的一组数据：,,,,27,,解：使 最小。

对应的,超定方程组为：,,,,即,,,,,,Matlab实现： u = A \ y ( 见 P228 例7.8 ),,,记为 Au = y,,28,,§,5.,数值积分,,复化梯形公式的实现:,,,,,trapz (y) * h,,即,,y = [ f(x,0,), … , f(x,n,) ],,29,,例: P233 例7.10,,%fun.m,,function 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) * d,,z =,,0.988,,,被积函数,,30,,§,5.,常微分方程数值解法,以Euler法为例例:,,,,,取步长 h = 0.2,,31,,% f.m,,function r = f (x, y),,r = y- x+ 1 ;,,,% orig.m,,function r = orig (x, y),,r = x+ exp (x) ;,,右端函数,准确解,,32,,% euler.m Euler方法的函数,,function [ X, Y ] = euler (a, b, h, c) % c= y,0,,n = round ((b-a)/h) + 1 ; % 计算点的个数,,x = zeros (n, 1) ; % 设定 x 维数,,y = zeros (n, 1) ; % 设定 y 维数,,x(1) = a ;,,y(1) = c ; % 初始条件,,str = sprintf (‘x0=%g, y0=%g \ n’, x(1), y(1));,,disp (str); % 显示初始条件,,33,,% 对每个点用,Euler公式计算,,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); % 显示计算结果与误差,,end,,X = x ;,,Y = y ; % 返回计算结果,,34,,% 主文件 main.m,,clear ; % 清除内存变量,,a=0; b=5; h=0.2; c=1; % 输入参数值,,[ X, Y ] = euler (a, b, h, c) ; %调用Euler函数,,% 绘图,,plot ( X, orig(X), ‘b-’, X, Y, ‘r--’) ;,,legend ( ‘精确解’, ‘数值解’) ;,,,,35,,>> main,,x0=0, y0=1,,1 x= 0.2 y= 1.4 e= 0.0214028,,2 x= 0.4 y= 1.84 e= 0.0518247,,……,,,25 x= 5 y= 100.396 e= 53.0169,,,,,,,,36,,。