matlab7教程课件第8章计算方法的matlab实现.ppt

MATLAB 7.0从入门到精通,主要讲述内容,第1章 MATLAB简介 第2章 数值运算 第3章 单元数组和结构 第4章 字符串 第5章 符号运算 第6章 MATLAB绘图基础 第7章 程序设计 第8章 计算方法的MATLAB实现 第9章 优化设计 第10章 Simulink仿真初探,第8章 计算方法的MATLAB实现,8.1 方程求根 roots见多项式求根；roots(多项式系数矩阵) fzero可求解非线性方程 fzero(function,x0)其中function为求解的方程，x0为估计的根，x0可为标量或长度为2的向量，为向量时函数的两端的值应该符号相反，此时求区间上的解。只能求解x0附近的一个解。即使在某个区间内有多个解，但是区间端点符号相同的话仍然出错。, fzero(x3-3*x-1,2) ans = 1.8794 fzero(x3-3*x-1,2,4) ? Error using = fzero The function values at the interval endpoints must differ in sign. fzero(x3-3*x-1,1,4) ans = 1.8794, fzero(x2-3*x+2,0) ans = 1.0000 fzero(x2-3*x+2,3) ans = 2.0000 fzero(x2-3*x+2,0,4) ? Error using = fzero The function values at the interval endpoints must differ in sign.,8.2 线性方程组数值解法 1、直接解法 A*X=B X=AB或x=inv(A)*B 例： 12x(1)-3x(2)+3x(3)=15 -18x(1)+3x(2)-x(3)=-15 X(1)+x(2)+x(3)=6, a=12 -3 3;-18 3 -1;1 1 1; b=15;-15;6; x=ab x = 1.0000 2.0000 3.0000 inv(a)*b ans = 1 2 3,2、线性方程组求解中的变换 上三角变换 U=triu(x)返回矩阵x的上三角部分； U=triu(x,k)返回第k条对角线以上部分的元素。, a=12 -3 3;-18 3 -1;1 1 1; triu(a) ans = 12 -3 3 0 3 -1 0 0 1, triu(a,1) ans = 0 -3 3 0 0 -1 0 0 0, triu(a,-1) ans = 12 -3 3 -18 3 -1 0 1 1,对角变换 U=diag(x)返回矩阵x主对角线上的元素，返回结果是一列向量形式； U=diag(x,k)返回第k条对角线上的元素值。 当x为向量时生成矩阵。, a=12 -3 3;-18 3 -1;1 1 1; diag(a) ans = 12 3 1, a=12 -3 3;-18 3 -1;1 1 1; diag(a,1) ans = -3 -1, a=12 -3 3;-18 3 -1;1 1 1; diag(a,-1) ans = -18 1,下三角变换 U=tril(x)返回矩阵x的下三角部分； U=tril(x,k)返回第k条对角线以上下部分的元素。, a=12 -3 3;-18 3 -1;1 1 1; tril(a) ans = 12 0 0 -18 3 0 1 1 1, tril(a,1) ans = 12 -3 0 -18 3 -1 1 1 1, tril(a,-1) ans = 0 0 0 -18 0 0 1 1 0,3、迭代解法 迭代解法非常适用于求解大型稀疏系数矩阵的方程组，在线性方程组常用的迭代解法主要有Jacobi迭代法、Gauss-Seidel迭代法。,Jacobi迭代法 编写jacobi函数，并使之计算 a=10 -2 -1;-2 10 -1;-1 -2 5; b=3 15 10; s=jacobi(a,b,0 0 0,eps) s = 1 2 3,function s=jacobi(a,b,x0,eps) %jacobi迭代法皆线性方程组 %a为系数矩阵，b为方程组ax=b中的右边的矩阵b，x0为迭代初值 if nargin=3 eps=1.0e-6; elseif nargin=eps x0=s; s=B*x0+f; end return,Gauss-Saidel迭代法 编写gauss函数，并使之计算 a=10 -2 -1;-2 10 -1;-1 -2 5; b=3 15 10; s=gauss(a,b,0 0 0,eps) s = 1 2 3,function s=gauss(a,b,x0,eps) %gauss-seidel迭代法皆线性方程组 %a为系数矩阵，b为方程组ax=b中的右边的矩阵b，x0为迭代初值 if nargin=3 eps=1.0e-6; elseif nargin=eps x0=s; s=B*x0+f; end return,8.3 非线性方程组数值解法 与线性方程组的求解一样，非线性方程组的求解也是应用很广泛的课题。一般情况，非线性方程组的数据值解法主要采用迭代法来求解。比较常用的迭代法主要有不动点迭代法、Newton迭代法、拟Newton迭代法等几种方法。,不动点迭代 编写不动点迭代函数，并使之计算 用不动点迭代法求解下面的方程组,首先编写上述非线性方程组的M文件fx.m staticiterate(0 0) ans = 1.0000 1.0000,function s=staticiterate(x,eps) %不动点迭代法解非线性方程组，x为迭代初值，eps为允许误差 if nargin=1 eps=1.0e-6; elseif nargin=eps x=xx; xx=fx(x); end s=xx; return,function y=fx(x) y(1)=0.1*(x(1)*x(1)+x(2)*x(2)+8); %y(1)=x(1)*x(1)+x(2)*x(2)-10*x(1)+8; y(2)=0.1*(x(1)*x(2)*x(2)+x(1)+8); %y(2)=x(1)*x(2)*x(2)+x(1)-10*x(2)+8; y=y(1) y(2);,Newton迭代法 用newton迭代法求解下面的方程组,首先，编写上述非线性方程组的M文件fx1.m 然后，编写上述非线性方程组导数的M文件dfx1.m newtoniterate(0 0) ans = 1 1,function s=newtoniterate(x,eps) %newton迭代法解非线性方程组，x为迭代初值，eps为允许误差 if nargin=1 eps=1.0e-6; elseif nargin=eps x=x0+x; x1=fx1(x); x2=-dfx1(x); x3=inv(x2); x0=x3*x1; end s=x0+x; return,function y=fx1(x) y(1)=x(1)*x(1)-10*x(1)+x(2)*x(2)+8; y(2)=x(1)*x(2)*x(2)-10*x(2)+x(1)+8; y=y(1) y(2);,function y=dfx1(x) y(1)=2*x(1)-10; y(2)=2*x(2); y(3)=x(2)*x(2)+1; y(4)=2*x(1)*x(2)-10; y=y(1) y(2);y(3) y(4);,8.4 插值与拟合 1、一维线性插值 yi=interp1(x,y,xi)返回在插值向量xi处的函数向量yi，它是根据向量x和y插值而来。如y是矩阵，则对y每一列进行插值，如xi中元素不在x内，返回NaN。 yi=interp1(y,xi)省略x，表示x=1:N，此时N为向量y的长度或为矩阵y的行数。 yi=interp1(x,y,xi,method)表示用method指定的插值方法进行插值。,Method可取如下的值： nearest最近插值 linear线性插值 spline三次样条插值 cubic三次插值 Method默认值为线性插值，上述插值要求向量x单调。, x=1 2 4 6 8 9 10 13 15 16; y=5 7 8 10 13 14 15 17 19 20; x1=1.2 2.1 3; y1=interp1(x,y,x1) y1 = 5.4000 7.0500 7.5000, x=1 2 4 6 8 9 10 13 15 16; y=5 7 8 10 13 14 15 17 19 20; x1=1.2 2.1 3; y1=interp1(x,y,x1,linear) y1 = 5.4000 7.0500 7.5000, x=1 2 4 6 8 9 10 13 15 16; y=5 7 8 10 13 14 15 17 19 20; x1=1.2 2.1 3; y1=interp1(x,y,x1,nearest) y1 = 5 7 8, x=1 2 4 6 8 9 10 13 15 16; y=5 7 8 10 13 14 15 17 19 20; x1=1.2 2.1 3; y1=interp1(x,y,x1,spline) y1 = 5.5520 7.1114 7.6785, x=1 2 4 6 8 9 10 13 15 16; y=5 7 8 10 13 14 15 17 19 20; x1=1.2 2.1 3; y1=interp1(x,y,x1,cubic) y1 = 5.5006 7.0814 7.5476,x=linspace(0,2*pi,100); y=sin(x); x1=0.5 1.4 2.1 2.6 4.2; y1=interp1(x,y,x1,nearest) plot(x,y,x1,y1,+) title(nearest插值) y1 = 0.4862 0.9848 0.8660 0.5137 -0.8660,x=linspace(0,2*pi,100); y=sin(x); x1=0.5 1.4 2.1 2.6 4.2; y1=interp1(x,y,x1,linear) plot(x,y,x1,y1,+) title(linear插值) y1 = 0.4793 0.9853 0.8631 0.5155 -0.8713,x=linspace(0,2*pi,100); y=sin(x); x1=0.5 1.4 2.1 2.6 4.2; y1=interp1(x,y,x1,spline) plot(x,y,x1,y1,+) title(spline插值) y1 = 0.4794 0.9854 0.8632 0.5155 -0.8716,x=linspace(0,2*pi,100); y=sin(x); x1=0.5 1.4 2.1 2.6 4.2; y1=interp1(x,y,x1,cubic) plot(x,y,x1,y1,+) title(cubic插值) y1 = 0.4794 0.9854 0.8632 0.5155 -0.8716,2、二维线性插值 zi=interp2(x,y,z,xi,yi)返回在插值向量xi、yi处的函数值向量，它是根据向量x、y与z插值而来，这里的x、y与z也可以是矩阵形式。如果xi、yi有元素不在x、y范围内，则返回NaN。 zi=interp2(z,xi,yi)省略x、y，表示x=1:N，y=1:M。此处M,N=size(z)。 zi=interp2(z,ntimes)表示在z的各点之间插入数据点来扩展z，依次执行ntimes次迭代，ntimes默认为1。 zi=interp2(x,y,z,xi,yi,method)表示用method指定的插值方法进行插值。Method取值同上，要求x与y都单调且具有相同格式。, z1=3 5 7 9 11 10 4 15;1 3 2 6 11 5 7 13; z=interp2(z1,1 3,1 2) z = 3 2, z1=3 5 7 9 11 10 4 15;1 3 2 6 11 5 7 13; z=interp2(z1)或z=interp2(z1,1) z = Columns 1 through 12 3.0000 4.0000 5.0000 6.0000 7.0000 8.0000 9.0000 10.0000 11.0000 10.5000 10.0000 7.0000 2.0000 3.0000 4.0000 4.2500 4.5000 6.0000 7.5000 9.2500 11.0000 9.2500 7.5000 6.5000 1.0000 2.0000 3.0000 2.5000 2.0000 4.0000 6.0000 8.5000 11.0000 8.0000 5.0000 6.0000 Columns 13 through 15 4.0000 9.5000 15.0000 5.5000 9.7500 14.0000 7.0000 10.0000 13.0000, x,y,z=peaks(10); meshc(x,y,z), x,y,z=peaks(10); xi,yi=meshgrid(-3:0.1:3,-3:0.1:3); zi=interp2(x,y,z,xi,yi,neareat); meshc(xi,yi,zi) title(nearest插值), x,y,z=peaks(10); xi,yi=meshgrid(-3:0.1:3,-3:0.1:3); zi=interp2(x,y,z,xi,yi,linear); meshc(xi,yi,zi) title(linear插值), x,y,z=peaks(10); xi,yi=meshgrid(-3:0.1:3,-3:0.1:3); zi=interp2(x,y,z,xi,yi,spline); meshc(xi,yi,zi) title(spline插值), x,y,z=peaks(10); xi,yi=meshgrid(-3:0.1:3,-3:0.1:3); zi=interp2(x,y,z,xi,yi,cubic); meshc(xi,yi,zi) title(cubic插值),x,y=meshgrid(-3:0.1:3,-3:0.1:3); z=peaks(x,y); meshc(x,y,z) title(原始图),3、三维线性插值 vi=interp3(x,y,z,v,xi,yi,zi)返回三维函数v在插值向量xi，yi，zi处的函数向量vi，它们大小形式相同。 vi=interp3(v,xi,yi,zi)省略x、y、z。同前。 vi=interp3(x,y,z,v,xi,yi,zi,method)同前。, x,y,z,v=flow(10); xi,yi,zi=meshgrid(0.1:0.25:10,-3:0.25:3,-3:0.25:3); vi=interp3(x,y,z,v,xi,yi,zi); slice(xi,yi,zi,vi,6 9.5,2,-2,0.2),4、Lagrange插值,编写拉氏插值函数然后运行 x=1 2 3 4 5; y=2 4 6 8 10; s=lagrange(x,y,1.6) s = 3.2000 s=lagrange(x,y,1.4 2.5 3.7) s = 0 0 7.4000 用法： s=lagrange(x,y,x0),5、newton插值,编写牛顿插值函数然后运行 x=0.4 0.55 0.65 0.8 0.9 1.05; y=0.41075 0.57815 0.69675 0.88811 1.02652 1.25382; s=newton(x,y,0.596,4) s = 0.6319 用法： s=newton(x,y,x0,n)其中n为插值多项式的次数。,6、三次样条插值 yy=spline(x,y,xx)利用三次样条插值法寻找在插值点xx处的插值函数 x=0:10; y=sin(x); xx=0:0.25:10; yy=spline(x,y,xx); plot(x,y,o,xx,yy),7、最小二乘法曲线拟合 p=polyfit(x,y,n)返回拟合的多项式的系数，系数存储在向量p中。 p,s=polyfit(x,y,n)返回是拟合生成的多项式系数向量p及用polyval函数获得的误差预测结果s。 x=1 3 4 5 6 7 8 9 10; y=10 5 4 2 1 1 2 3 4 ; p,s=polyfit(x,y,4); y1=polyval(p,x); plot(x,y,go,x,y1,b-),8.5 数值积分 1、Newton-Cotes求积公式 梯形积分法 s=trapz(x,y)表示用梯形法求y对x的积分。 x=0:0.1:pi; y=sin(x); s=trapz(x,y) s = 1.9975, x=0:0.01:pi; y=sin(x); s=trapz(x,y) s = 2.0000 由此可见，积分精度与步长有关。,simpson法数值积分 此方法的数值积分用函数quad来实现。 q=quad(f,a,b)表示使用自适应递归的simpson方法从区间a到b对函数f(x)积分。求积相对误差在1e-3范围内。 q=quad(f,a,b,tol)表示使用自适应递归的simpson方法从区间a到b对函数f(x)积分。求积相对误差在tol范围内。 此函数最大递归调用次数为10次，如超出则返回Inf。, quad(sin(x),0,pi) ans = 2.0000,cotes(科特斯)法数值积分 此方法的数值积分用函数quad8来实现。 q=quad8(f,a,b)表示使用自适应递归的newton-cotes8方法从区间a到b对函数f(x)积分。求积相对误差在1e-3范围内。 q=quad8(f,a,b,tol)表示使用自适应递归的newton-cotes8方法从区间a到b对函数f(x)积分。求积相对误差在tol范围内。 此函数最大递归调用次数为10次，如超出则返回Inf。, quad(cos(x),0,pi/2) ans = 1.0000,2、Gauss求积公式 首先编写4点gauss求积公式的M文件，然后运行。 用编写的函数求积分 首先编写求积函数gaussff.m, gausslegendre(0,1) ans = 0.7183,3、romberg求积公式 编写romberg积分函数。求下函数的积分： 首先编写积分函数rombergff.m, romberg(0,1) ans = 0.6366,4、二重积分 q=dblquad(fun,x1,x2,y1,y2) q=dblquad(fun,x1,x2,y1,y2,tol)tol指定精度 q=dblquad(fun,x1,x2,y1,y2,tol,method) q=dblquad(3*x.2+3*y.2,0,1,0,1) q = 2,5、三重积分 q=triplequad(fun,x1,x2,y1,y2,z1,z2) q=triplequad(fun,x1,x2,y1,y2,z1,z2,tol)tol指定精度 q=triplequad(fun,x1,x2,y1,y2,z1,z2,tol,method) q=triplequad(x.2+y.2+z.2,0,1,0,1,0,1) q = 1,8.6 常微分方程的数值解法 1、Euler方法 编写euler法的M函数，编写搜求函数文件eulerff.m，然后求解 y=-2xy 0 euler(eulerff,0,1.2,1,100) ans = 0.2370 用法： euler(fun,x0,xn,y0,n)n为区间分成的份数。,2、Runge-kutta法