数值分析：Ch10.数值算法实现.ppt

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,Ch10.,数值算法实现,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,、,顺序,Gauss,消去法,（,1,）,消去过程：,第,k,步，计算,（,2,）回代过程：,%,文件,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),1e6,a(i,j)=a(i,j),mult,*a(k,j);,%else,%,disp,(,顺序,Gauss,消去法失败,);,%,pause,；,%exit,；,%end,end,可,去掉%，,判断主元是否为0,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;,回代,例：,%,主文件,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,3,、,Jacobi,迭代法,%,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,例,:,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,解,向量,迭代次数,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),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=,56,ans,=,1.6180,先,编写函数,2,、不动点,迭代法,%,iterate.m,function y=iterate(x),x1=g(x);,n=1;,while (abs(x1 x)=1.0e 6)&(n iterate(3),x1=,3.7331,n=,22,初值不同,迭代步数会不同,迭代函数,3,、牛顿,迭代法,%,newton.m,function 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.7331,n=n=,6 12,比,不动点,迭代快,习题三,1,、编写下三角形方程组前推法的函数；,2,、编写,Jacobi,迭代分量形式,的函数；,3,、利用牛顿,迭代法,求解,计算方法,P114,习题,4.7,，写出各函数和命令文件。,下 次 见,3.,插值方法,以牛顿插值为例。,已知函数 在互异点 上的值为,。,n,次,Newton,插值多项式,为：,定义,设函数 在互异点 上的值为,。定义：,1,）,一阶均差,为,2,）,二阶均差,为,3,）递推下去，,k,阶均差,为,例：,已知,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,处的,值。比较与原函数的误差。,%,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 阶均差,%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 阶均差,%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);,插值函数,%,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(,原函数,插值函数,);,main,x=0.75,处值与误差,l=,0.9221,e=,0.0150,x=1.75,处值与误差,l=,2.4228,e=,0.0077,4.,曲线拟合的最小二乘法,1,、多项式拟合,2,、一般情形,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,多项式次数,例,:,用三次多项式在,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 ),拟合函数,2,、一般情形,对应于解超定方程组,用矩阵除法实现,例：用,y=a+bx,2,拟合给定的一组数据：,解：使 最小。,对应的,超定方程组为：,即,Matlab,实现：,u=A y (,见,P228,例,7.8),记为,Au=y,5.,数值积分,复化梯形公式的实现,:,trapz,(y)*h,即,y=f(x,0,),f(,x,n,),例,:,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.90084027660688,被积,函数,5.,常微分方程数值解法,以,Euler,法为例。,例,:,取步长,h=0.2,%,f.m,function r=f(x,y),r=y-x+1;,%orig.m,function r=,orig,(x,y),r=x+exp(x);,右端函数,准确解,%,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,);%,显示初始条件,%,对每个点用,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;%,返回计算结果,%,主文件,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(,精确解,数值解,);,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,谢 谢,See You Next Time!,