计算方法书中的MATLAB实例PPT学习课件.ppt

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,#,用,Matlab,实现求积分,在,Matlab,程序编辑器中输入：,y0=log(6.0/5.0);,fprintf(y%d=%fn,0,y0),n=1;,while(1),yl=1.0/n-5*y0;,fprintf(y%d=%fn,n,yl),if(n=20)break;,end,y0=yl;,n=n+1;,end,1,2025/3/25 周二,2,2025/3/25 周二,程序运行结果,:,3,2025/3/25 周二,function x=nabisect(fname,a,b,e)%fname,为内嵌函数,表达式；,a,b,为区间端点；,e,为输入定义的精度,if nargin0,error(,函数在两端点值必须异号,);end,x=(a+b)/2,while(b-a)(2*e),fx=feval(fname,x);,if fa*fx0,b=x;fb=fx;,else a=x;fa=fx;end,x=(a+b)/2,end,程序如下：,Matlab,应用举例,(,二分法,),4,2025/3/25 周二,5,2025/3/25 周二,程序运行结果：,6,2025/3/25 周二,用,Matlab,实现,Newton,迭代法求解,function x=nanewton(fname,dfname,x0,e,N),if nargin5,N=500;end%N,为迭代次数,if nargine&kN,k=k+1;x0=x;,x=x0-feval(fname,x0)/feval(dfname,x0);,disp(x),end,if k=N,warning(,已达迭代次数上限,);end,程序如下,：,7,8,2025/3/25 周二,程序运行结果：,9,2025/3/25 周二,用,Matlab,实现顺序,Gauss,消去法,在,Matlab,程序编辑器中输入：,function x=nagauss(a,b,flag)%,解线形方程组,ax=b,，,a,为系数矩阵，,b,为右端列向量，,flag,若为,0,，则显示中间过程，否则不显示，默认为,0,，,x,为解向量,if nargin3,flag=0;end,n=length(b);a=a,b;,%,消元,for k=1:(n-1),a(k+1):n,(k+1):(n+1)=a(k+1):n,(k+1):(n+1)-a(k+1):n,k)/a(k,k)*a(k,(k+1):(n+1);,a(k+1):n,k)=zeros(n-k,1);,if flag=0,a,end,end,10,2025/3/25 周二,%,回代,x=zeros(n,1);,x(n)=a(n,n+1)/a(n,n);,for k=n-1:-1:1,x(k,:)=(a(k,n+1)-a(k,(k+1):n)*x(k+1):n)/a(k,k);,end,11,2025/3/25 周二,12,2025/3/25 周二,程序运行,结果：,13,2025/3/25 周二,用,Matlab,实现选列主元,Gauss,消去法解线性方程组,在,Matlab,程序编辑器中输入：,function x=nagauss2(a,b,flag)%a,为系数矩阵；,b,为右端列向量；,flag,若为,0,，则显示中间过程，否则不显示,if nargink,t=a(k,:);a(k,:)=a(p,:);a(p,:)=t;end,%,消元,a(k+1):n,(k+1):(n+1)=a(k+1):n,(k+1):(n+1)-a(k+1):n,k)/a(k,k)*a(k,(k+1):(n+1);,a(k+1):n,k)=zeros(n-k,1);,14,2025/3/25 周二,if flag=0,a,end,end,%,回代,x=zeros(n,1);,x(n)=a(n,n+1)/a(n,n);,for k=n-1:-1:1,x(k,:)=(a(k,n+1)-a(k,(k+1):n)*x(k+1):n)/a(k,k);,end,15,2025/3/25 周二,16,2025/3/25 周二,程序运行,结果,:,17,2025/3/25 周二,用,Matlab,实现,LU,分解,在,Matlab,程序编辑器中输入：,function L,U=nalu(a)%a,为可逆方阵；,L,返回单位下三角矩阵；,U,返回上三角矩阵,n=length(a);,U=zeros(n,n);L=eye(n,n);,U(1,:)=a(1,:);L(2:n,1)=a(2:n,1)/U(1,1);,for k=2:n,U(k,k:n)=a(k,k:n)-L(k,1:k-1)*U(1:k-1,k:n);,L(k+1:n,k)=(a(k+1:n,k)-L(k+1:n,1:k-1)*U(1:k-1,k)/U(k,k);,end,18,2025/3/25 周二,19,2025/3/25 周二,程序运行,结果,:,20,2025/3/25 周二,用,Matlab,实现,Lagrange,插值：,function yy=nalagr(x,y,xx)%x,为节点；,y,为节点值；,xx,为插值点；,yy,为返回值,m=length(x);n=length(y);,if m=n,error(,向量,x,与,y,长度必须一致,);end,s=0;,for i=1:n,t=ones(1,length(xx);,for j=1:n,if j=i,t=t.*(xx-x(j)/(x(i)-x(j);,end,end,s=s+t*y(i);,end,yy=s;,21,2025/3/25 周二,22,2025/3/25 周二,求解例题得,一次插值和,二次插值,如右图所示：,23,2025/3/25 周二,用,Matlab,实现最小二乘拟合,function p=nafit(x,y,m)%x,y,为数据；,m,为拟合多项式次数,A=zeros(m+1,m+1);%zeros(),生成零矩阵,for i=0:m,for j=0:m,A(i+1,j+1)=sum(x.(i+j);,end,b(i+1)=sum(x.i.*y);,end,a=Ab;,p=fliplr(a);%fliplr(),为矩阵左右翻转,在,Matlab,程序编辑器中输入：,24,2025/3/25 周二,25,2025/3/25 周二,最小二乘,拟合多项式,结果如右图,:,26,2025/3/25 周二,用,Matlab,实现用,Euler,格式解常微分方程,在,Matlab,程序编辑器中输入：,function x,y=naeuler(dyfun,xspan,y0,h)%dyfun,为导函数；,xspan,为求解区间；,y0,为初值；,h,为步长,x=xspan(1):h:xspan(2);,y(1)=y0;,for n=1:length(x)-1,y(n+1)=y(n)+h*feval(dyfun,x(n),y(n);,end,x=x;y=y;,27,2025/3/25 周二,28,2025/3/25 周二,Euler,格式,求解得：,29,2025/3/25 周二,用,Matlab,实现用改进,Euler,格式解常微分方程,在,Matlab,程序编辑器中输入：,function x,y=naeuler2(dyfun,xspan,y0,h)%dyfun,为导函数；,xspan,为求解区间；,y0,为初值；,h,为步长,x=xspan(1):h:xspan(2);y(1)=y0;,for n=1:length(x)-1,k1=feval(dyfun,x(n),y(n);,y(n+1)=y(n)+h*k1;,k2=feval(dyfun,x(n+1),y(n+1);,y(n+1)=y(n)+h*(k1+k2)/2;,end,x=x;y=y;,30,2025/3/25 周二,31,2025/3/25 周二,用改进,Euler,格式求得：,32,2025/3/25 周二,