计算方法书中的MATLAB实例

1、用Matlab实现求积分 dx x x I n n 1 0 5 在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 计算方法书中的MATLAB实例 计算方法书中的MATLAB实例 程序运行结果: 计算方法书中的MATLAB实例 function x=nabisect(fname,a,b,e) %fname为内嵌函数 表达式；a,b为区间端点；e为输入定义的精度

2、 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实例 计算方法书中的MATLAB实例 程序运行结果： 计算方法书中的MATLAB实例 用Matlab实现Newton迭代法求解 位有效数字。内的一个实根，精确至，在区间 用迭代法求方程例 410 01)1( 5 2 xxf(x) function x=nanewton(fname,dfnam

3、e,x0,e,N) if nargin5,N=500;end % N为迭代次数 if narginex0=x; x=x0-feval(fname,x0)/feval(dfname,x0); disp(x) end if k=N,warning(已达迭代次数上限);end 程序如下： 计算方法书中的MATLAB实例 计算方法书中的MATLAB实例 程序运行结果： 计算方法书中的MATLAB实例 用用MatlabMatlab实现顺序实现顺序GaussGauss消去法消去法 在在MatlabMatlab程序编辑器中输入：程序编辑器中输入： function x=nagauss(a,b,flag) %

4、function x=nagauss(a,b,flag) %解线形方程组解线形方程组ax=bax=b， a a为系数矩阵，为系数矩阵，b b为右端列向量，为右端列向量，flagflag若为若为0 0，则显示中间，则显示中间 过程，否则不显示，默认为过程，否则不显示，默认为0 0，x x为解向量为解向量 if nargin3,flag=0;endif nargin3,flag=0;end n=length(b); a=a,b;n=length(b); a=a,b; % % 消元消元 for k=1:(n-1)for k=1:(n-1) a(k+1):n,(k+1):(n+1)=a(k+1):n,

5、(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+1):n,k)/a(k,k)* *a(k,(k+1):(n+1);a(k,(k+1):(n+1); a(k+1):n,k)=zeros(n-k,1); a(k+1):n,k)=zeros(n-k,1); if flag=0,a,end if flag=0,a,end end end 计算方法书中的MATLAB实例 % % 回代回代 x=zeros(n,1);x=zeros(n,1); x(n)=a(n,n+1)/a(n,n);x(n)=a(n

6、,n+1)/a(n,n); for k=n-1:-1:1for k=n-1:-1:1 x(k,:)=(a(k,n+1)- x(k,:)=(a(k,n+1)- a(k,(k+1):n)a(k,(k+1):n)* *x(k+1):n)/a(k,k);x(k+1):n)/a(k,k); endend 计算方法书中的MATLAB实例 计算方法书中的MATLAB实例 程序运行程序运行 结果：结果： 计算方法书中的MATLAB实例 用用MatlabMatlab实现选列主元实现选列主元GaussGauss消去法解线性方程组消去法解线性方程组 在在MatlabMatlab程序编辑器中输入：程序编辑器中输入：

7、function x=nagauss2(a,b,flag) %afunction x=nagauss2(a,b,flag) %a为系数矩阵；为系数矩阵；b b为为 右端列向量；右端列向量；flagflag若为若为0 0，则显示中间过程，否则不显示，则显示中间过程，否则不显示 if nargin3,flag=0;endif nargink,t=a(k,:); a(k,:)=a(p,:); a(p,:)=t; endif pk,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

8、+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+1):n,k)/a(k,k)* *a(k,(k+1):(n+1);a(k,(k+1):(n+1); a(k+1):n,k)=zeros(n-k,1); a(k+1):n,k)=zeros(n-k,1); 计算方法书中的MATLAB实例 if flag=0,a,endif flag=0,a,end end end % %回代回代 x=zeros(n,1);x=zeros(n,1); x(n)=a(n,n+1)/a(n,n);x(n)=a(n,n+1)/a(n,

9、n); for k=n-1:-1:1for k=n-1:-1:1 x(k,:)=(a(k,n+1)- x(k,:)=(a(k,n+1)- a(k,(k+1):n)a(k,(k+1):n)* *x(k+1):n)/a(k,k);x(k+1):n)/a(k,k); endend 计算方法书中的MATLAB实例 计算方法书中的MATLAB实例 程序运行程序运行 结果结果: : 计算方法书中的MATLAB实例 用用MatlabMatlab实现实现LULU分解分解 在在MatlabMatlab程序编辑器中输入：程序编辑器中输入： function L,U=nalu(a) % afunction L,U=

10、nalu(a) % a为可逆方阵；为可逆方阵；L L返回单位下返回单位下 三角矩阵；三角矩阵；U U返回上三角矩阵返回上三角矩阵 n=length(a);n=length(a); U=zeros(n,n);L=eye(n,n);U=zeros(n,n);L=eye(n,n); U(1,:)=a(1,:);L(2:n,1)=a(2:n,1)/U(1,1);U(1,:)=a(1,:);L(2:n,1)=a(2:n,1)/U(1,1); for k=2:nfor k=2:n U(k,k:n)=a(k,k:n)-L(k,1:k-1) U(k,k:n)=a(k,k:n)-L(k,1:k-1)* *U(1

11、:k-1,k:n);U(1:k-1,k:n); L(k+1:n,k)=(a(k+1:n,k)-L(k+1:n,1:k-1) L(k+1:n,k)=(a(k+1:n,k)-L(k+1:n,1:k-1)* *U(1:k-U(1:k- 1,k)/U(k,k);1,k)/U(k,k); endend 计算方法书中的MATLAB实例 计算方法书中的MATLAB实例 程序运行程序运行 结果结果: : 计算方法书中的MATLAB实例 用用MatlabMatlab实现实现 LagrangeLagrange插值：插值： function yy=nalagr(x,y,xx) % xfunction yy=nala

12、gr(x,y,xx) % x为节点；为节点；y y为节点值；为节点值； xxxx为插值点；为插值点；yyyy为返回值为返回值 m=length(x);n=length(y);m=length(x);n=length(y); if m=n,error(if m=n,error(向量向量x x与与y y长度必须一致长度必须一致);end);end s=0;s=0; for i=1:nfor i=1:n t=ones(1,length(xx); t=ones(1,length(xx); for j=1:n for j=1:n if j=i, if j=i, t=t. t=t.* *(xx-x(j)/

13、(x(i)-x(j);(xx-x(j)/(x(i)-x(j); end end end end s=s+t s=s+t* *y(i);y(i); endend yy=s;yy=s; 计算方法书中的MATLAB实例 计算方法书中的MATLAB实例 求解例题得求解例题得 一次插值和一次插值和 二次插值二次插值 如右图所示：如右图所示： 计算方法书中的MATLAB实例 用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+

14、1,j+1)=sum(x.(i+j); end b(i+1)=sum(x.i.*y); end a=Ab; p=fliplr(a); % fliplr()为矩阵左右翻转 在Matlab程序编辑器中输入： 计算方法书中的MATLAB实例 计算方法书中的MATLAB实例 最小二乘 拟合多项式 结果如右图: 计算方法书中的MATLAB实例 用用MatlabMatlab实现用实现用EulerEuler格式解常微分方程格式解常微分方程 在在MatlabMatlab程序编辑器中输入：程序编辑器中输入： function x,y=naeuler(dyfun,xspan,y0,h) %dyfunfunctio

15、n x,y=naeuler(dyfun,xspan,y0,h) %dyfun为为 导函数；导函数；xspanxspan为求解区间；为求解区间；y0y0为初值；为初值；h h为步长为步长 x=xspan(1):h:xspan(2);x=xspan(1):h:xspan(2); y(1)=y0;y(1)=y0; for n=1:length(x)-1for n=1:length(x)-1 y(n+1)=y(n)+h y(n+1)=y(n)+h* *feval(dyfun,x(n),y(n);feval(dyfun,x(n),y(n); endend x=x;y=y;x=x;y=y; 计算方法书中的

16、MATLAB实例 计算方法书中的MATLAB实例 EulerEuler格式格式 求解得：求解得： 计算方法书中的MATLAB实例 用用MatlabMatlab实现用改进实现用改进EulerEuler格式解常微分方程格式解常微分方程 在在MatlabMatlab程序编辑器中输入：程序编辑器中输入： function x,y=naeuler2(dyfun,xspan,y0,h) %dyfunfunction x,y=naeuler2(dyfun,xspan,y0,h) %dyfun 为导函数；为导函数；xspanxspan为求解区间；为求解区间；y0y0为初值；为初值；h h为步长为步长 x=xspan(1):h:xspan(2);y(1)=y0;x=xspan(1):h:xspan(2);y(1)=y0; for n=1:length(x)-1for n=1:length(x)-1 k1=feval(dyfun,x(n),y(n); k1=feval(dyfun,x(n),y(n); y(n+1)=y(n)+h y(n+1)=y(n)+h* *k1;k1; k2=feval(dyfun,x(