Matlab数值算法实现

1、Ch10. 数值算法实现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）消去过程：第 k 步，计算（2）回代过程： 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、数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先编写函数122、不动点迭代法 % 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初值不同,迭代步数会不同迭代函数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、比不动点迭代快163. 插值方法 以牛顿插值为例。 已知函数 在互异点 上的值为 。 n次Newton插值多项式 为：17定义 设函数 在互异点 上的值为 。定义：1）一阶均差为2）二阶均差为3）递推下去，k 阶均差为18例： 已知 f (x) = sqrt (2x) + ln (x) 节点 x1 = 0.5 , x2 = 1 , x3 = 1.5 , x4 = 2.0 求 Newton 插值多项式，并计算其在 x = 0.75 , x = 1.75处的值。比较与原函数的误差。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) ;原函数0 阶均差 (原函数)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) ;插值函数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 (原函数, 插值函数);23 mainx=0.75处值与误差l = 0.9221e = 0.0150 x=1.75处值与误差l = 2.4228e = 0.0077241、 多项式拟合 已知 ( xi , yi ) ( i = 1, , m) polyfit(x, y, n)x= x1 , x2 , , xm y= y1 , y2 , , ym 多项式次数25例: 用三次多项式在 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 )拟合函数262、一般情形对应于解超定方程组 用矩阵除法实现例：用 y = a + bx2 拟合给定的一组数据： 27解：使 最小。对应的超定方程组为：即Matlab实现： u = A y ( 见 P228 例7.8 ) 记为 Au = y285. 数值积分 复化梯形公式的实现: 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被积函数305. 常微分方程数值解法 以Euler法为例。例:取步长 h = 0.231% f.mfunction r = f (x, y)r = y- x+ 1 ;% orig.mfunction r = orig (x, y)r = x+ exp (x) ;右端函数准确解32% euler.m

13、Euler方法的函数function X, Y = euler (a, b, h, c) % c= y0n = 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); % 显示计算结果与误差endX = x ; Y = y ; % 返回计算结果34% 主文件 main.mclear ; % 清除内存变量a=0; b=5; h=0.2; c