Matlab数值计算.ppt

MATLAB与数值计算 刘东毅 数学与应用数学系 天津大学理学院 *1天津大学数学系 MATLAB在数值分析中的应用举例 线性代数方程组的数值解法 Gauss消取法解方程组； 多项式和矩阵的特征系统(eig) 函数的插值 Lagrange多项式插值 曲线拟合与函数的数值逼近 构造Legendre正交多项式 Date2天津大学数学系 MATLAB在数值分析中的应用举例 数值积分与数值微分 解算子quad, quadl。 计算椭圆积分 常微分方程（组）数值解 解初值问题的解算子 ODE23,ODE45,ODE113 非线性方程和方程组的数值解法 解算子FZERO，FSOLVE Date3天津大学数学系 【例A1.1】求解方程组 1.线性代数方程组的数值解法 这里 ， . Date4天津大学数学系 线性代数方程组的数值解法 (1) 在键盘上输入下列内容 A = 1,2,3; 4,5,6; 7,8,0%节尾没有分号； b = 37;85;61; %节尾有分号； x=Ab %节尾没有分号； (2) 每按一次【Enter】键，指令就被马上 执行（逐行执行）。由于第二条指令节 尾有分号，其结果不被显示出来，其它 两条指令的结果被马上显示出来。最后 在指令窗中将显示以下结果： Date5天津大学数学系 线性代数方程组的数值解法 A = x = 1 2 3 3.0000 4 5 6 5.0000 7 8 0 8.0000 B ACK Date6天津大学数学系 2.多项式和矩阵的特征系统 2.1多项式 MATLAB约定：用系数行向量 P=a0,a1,an-1,an 来表示多项式 如 P = 2 0 1 1 代表 Date7天津大学数学系 多项式的四则运算 多项式的加减法与一维数组的加减法类 似，只不过要注意多项式的阶数与行向 量元素个数的关系。 Date8天津大学数学系 乘除法：MATLAB提供了卷积和解卷函数 乘法：p = conv(p1,p2)，它表示多项式p 为多项式p1与多项式p2的积。 除法（带余除法）： q, r = deconv(p1, p2)，它表示多项式p1 被p2除的商为多项式q而余项是多项式r， 即满足p1 = q*p2+r。 Date9天津大学数学系 【例A.2.1】求下列多项式的“商”及“余” p1=conv(1,0,2,conv(1,4,1,1); p2=1 0 1 1; q,r=deconv(p1,p2); cq=商多项式为 ; cr=余多项式为 ; disp(cq,poly2str(q,s); disp(cr,poly2str(r,s); Date10天津大学数学系 2.2 多项式求根及其逆问题 n次多项式具有n个根（实根或成对的 共轭复根）。MATLAB提供的roots函数用 于求多项式的全部根，其调用格式为： x=roots(P) 其中P为多项式的系数向量，根赋给向量x ，即x(1),x(2),x(n)分别代表多项式的n个 根。反之，P = poly( x ) 函数poly生成以向量x为根的多项式. Date11天津大学数学系 【例A.2.2】多项式求根及其逆问题 R=1,-2, -0.3+0.5*i,-0.3-0.5*i; P=poly(R) PR=real(P) PPR=poly2str(PR,x) Rr=roots(P) Date12天津大学数学系 结果依次为 P = 1.0000 1.6000 -1.0600 -0.8600 -0.6800 PR = 1.0000 1.6000 -1.0600 -0.8600 -0.6800 PPR=x4+1.6 x3 - 1.06 x2 - 0.86 x - 0.68 Rr = -2.0000 1.0000 -0.3000 + 0.5000i -0.3000 - 0.5000i Date13天津大学数学系 2.3 矩阵的特征系统 【例A.2.3】计算的矩阵A全部特征值， 特征向量与特征多项式。 Date14天津大学数学系 矩阵的特征系统 解：在指令窗输入： A=3,-4,3;-4,6,3;3,3,1 V,D=eig(A) P = poly(A) 便可得到结果： V = -0.5818 0.6312 0.5130 -0.4534 0.2719 -0.8488 0.6752 0.7264 -0.1280 每一列为A的一个特征向量。 Date15天津大学数学系 矩阵的特征系统 D = -3.5995 0 0 0 4.7296 0 0 0 8.8699 D的对角线上的元素为对应的特征值，即 对应V的每一个列向量。 P =1.0000 -10.0000 -7.0000 151.0000 BACK Date16天津大学数学系 3.函数的插值 Lagrange多项式插值 【例A.3.1】已知 的三个数据点 ， 和 。求二次 Lagrange插值函数L2。 Date17天津大学数学系 函数的插值(c1) 解：编写计算L2的程序： clear all x=-1,0,1; y=4,-1,2; p,pexpr=lagrpoly(x,y); 在指令窗输入pexpr，可得到L2的 表达式：4 x2 1 x 1。 BACK Date18天津大学数学系 4.曲线拟合与函数的数值逼近 【例A.4.1】编写MATLAB程序，构造 上的Legendre正交多项式。 解：本题的目的是为用正交多项式进行函 数的最佳平方逼近提供正交多项式子函数 。同时进一步增加对正交多项式的感性认 识。该程序利用两个cell数组：正交多项 式的行向量形式（P）和正交多项式的字 符串形式（Pexpr）来存放多项式，程序 如下： Date19天津大学数学系 曲线拟合与函数的数值逼近(c1) P,Pexpr=plegendre(7,1); format rat celldisp(P); celldisp(Pexpr); format short g Date20天津大学数学系 曲线拟合与函数的数值逼近(c2) Date21天津大学数学系 曲线拟合与函数的数值逼近(c3) Date22天津大学数学系 曲线拟合与函数的数值逼近(c4) 可得出这八个正交多项式为：ep1 = 1； ep2 = x；ep3 = 1.5 x2 - 0.5； ep4 = 2.5 x3 - 1.5 x；ep5 = 4.375 x4 - 3.75 x2 + 0.375； ep6 = 7.875 x5 - 8.75 x3 + 1.875 x； ep7 = 14.4375 x6 - 19.6875 x4 + 6.5625 x2 - 0.3125； ep8 = 26.8125 x7 - 43.3125 x5 + 19.6875 x3 - 2.1875 x BACK Date23天津大学数学系 【例A.4.2】已知一组实验数据如下: 用二阶多项式曲线进行拟合。编写 MATLAB程序，计算此二阶多项式表达式 .并绘出拟合图形。 xk 00.250.500.751.00 yk 1.00001.28401.64872.11702.7183 Date24天津大学数学系 程序如下： %准备数据 x=0 0.25 0.50 0.75 1.00; y=1.0000,1.2840,1.6487,2.1170, 2.7183; %调用拟合函数 pp=polyfit(x,y,2); %将多项式行向量形式转成字符串形式 p2=poly2str(pp,x) Date25天津大学数学系 %绘图，折线图，准备数据点对(x,y) xi=0:0.01:1; yi=polyval(pp,xi);%计算多项式在xi处的值 %调用绘图函数：原始数据,二阶曲线 plot(x,y,o,xi,yi); %图形修饰 legend(原始数据,二阶曲线); title(多项式拟合) axis equal; axis equal; xlabel(x);ylabel(y=p(x); Date26天津大学数学系 Date27天津大学数学系 5.数值积分与数值微分 【例A.5.1】计算椭圆积分,控制精度106 。 解：程序如下： fun = inline(sqrt(1-0.8*sin(x).2) ,x); v1, fcn1 = quad(fun, 0, 2*pi); v2, fcn2 = quadl(fun, 0, 2*pi); Date28天津大学数学系 数值积分与数值微分(c1) quad 和quadl的精度控制输入宗量在第四个位 置，却省值为106。故在上述程序中 “ = quad(fun, 0, 2*pi)”等价于 “ = quad(fun, 0, 2*pi, 1e-6)”。 结果为： v14.7139593，被积函数估值次数fcn1 = 81； v24.7139597，被积函数估值次数fcn2 = 168。 BACK Date29天津大学数学系 6.常微分方程（组）数值解 MATLAB求解的一阶常微分方程（组 ）应具有形式(初值问题)： Date30天津大学数学系 或写成向量形式： 其中 Date31天津大学数学系 6.常微分方程（组）数值解 MATLAB为解决常微分方程初值问题 提供了一套系统的数值解法程序,它包括解 算子(Solver)、ODE文件和算法参数选项 options . 例如ODE45 ，ODE23 ， ODE15S， ODE113，ODE23S，ODE23T和ODESET 等，其它函数的使用请参见在线帮助或 MATLAB使用手册。 Date32天津大学数学系 6.常微分方程（组）数值解 解算子是指MATLAB提供的各种常微分 方程初值问题数值解法程序,如ode45和 ode15s等； ODE文件是指被解算子调用的，由用户 自己编写的，计算导数的函数f(t,y)的M-函 数文件（f(t,y) 也称为ODE函数）； Options 选项是可用odeset指令来设置一 些可选的参数值. Date33天津大学数学系 我们以ode45为例介绍解算子使用方法. 函数ODE45 用于求解非刚性系统其 调用格式为： T,Y = ODE45(F,TSPAN,Y0)， T,Y = ODE45(F,TSPAN,Y0,OPTIONS)， T,Y=ODE45(F,TSPAN,Y0,OPTIONS,P1,P2,. ) Date34天津大学数学系 说明： 这里TSPAN = T0,TF 微分系统 dy/dt = F(t,y) 的积分区间，Y0为初始条 件，F是一ODE函数文件名，此文件定义 了函数F(T,Y)，其必须返回一列向量 解矩阵Y的每一行对应于返回时间列 向量T的一元素要想得到在固定时间序 列T0, T1, ., TF（单增或单减）的解，使 用格式TSPAN = T0,T1, . ,TF。 Date35天津大学数学系 说明： options用来设置一些参数值，缺省时，相 当于调用格式1。如数量的相对误差 RelTol 缺 省值为1e-3，向量各分量缺省的绝对误差 AbsTol缺省值为1e-6具体参见odeset P1,P2,.的作用是传递附加参数P1,P2,.到 ODE文件，如：F(T,Y,FLAG,P1,P2,.)（参见 ODEFILE)；若没有OPTIONS参数设置，须在 OPTIONS参数设置位置上设置OPTIONS = ， 以便正确传递参数 Date36天津大学数学系 【例A.6.1】解如下常微分方程初值问题： 解：首先将方程（组）写成标准形式： Date37天津大学数学系 第二，编写描述问题的ODE函数文件 ： F(x,Y) function df = odefun(x,y) df = 3.*y./(1+x); 第三，调用解算子ODE45求解 tspan = 0,1;Y0 = 1; T4,Y4 = ode45(odefun,tspan,Y0); 算法采用变步长自适应的Runge-Kutta法 Date38天津大学数学系 Date39天津大学数学系 利用ODE23，ODE45和ODE113重新求解 ，以比较它们的区别。程序如下： clear all fun = inline(3.*y./(1+x); ftr = inline(1+x).3); tspan = 0,1;Y0 = 1; T2,Y2 = ode23(fun,tspan,Y0); T4,Y4 = ode45(fun,tspan,Y0); Ta,Ya = ode113(fun,tspan,Y0); err2 = Y2-feval(ftr,T2) err4 = Y4-feval(ftr,T4) erra = Ya-feval(ftr,Ta) Date40天津大学数学系 常微分方程（组）数值解(c2) subplot(221) fplot(ftr,0,1); hold on; plot(T2,Y2,.); hold off legend(true,ode23,4); axis(0,1,0,10) subplot(222) fplot(ftr,0,1) hold on; plot(T4,Y4,.); hold off; legend(true,ode45,4) Date41天津大学数学系 常微分方程（组）数值解(c3) subplot(223) fplot(ftr,0,1) hold on; plot(Ta,Ya,.); hold off; legend(true,ode113,4) subplot(224) plot(T2,err2,o,T4,err4,.,Ta,erra,d) legend(ode23,ode45,ode113,3); title(误差的比较) Date42天津大学数学系 常微分方程（组）数值解(c4) BACK Date43天津大学数学系 【例A.6.2】单摆系统模型。 运动方程为: 其中表示摆杆与垂直方向的夹角，单位 是弧度, c=m/k, m为摆锤的质量，k为摩擦 阻力系数，L为摆杆长度（假设摆杆质量 为零），g为重力加速度。设初始条件为 Date44天津大学数学系 L=0.5米，g=9.8米/秒，c=-0.05，=1 解：设 则 以上方程组即为标准形式。 Date45天津大学数学系 编写描述问题的ODE函数文件 ：F(x,Y) function y = odefunp(t,y) g = 9.8; L = 0.5; c = 0.1; y=y(2);-g./L.*sin(y(1)-c.*y(2); Date46天津大学数学系 编写调用解算子ODE45求解的主文件 function pendulum Y0=pi/4; 0 ; TSPAN=0; 20; T,Y=ode45(odefunp,TSPAN,Y0) plot(T,Y(:,1); figure plot(Y(:,1),Y(:,2); %相图 pendulum.m Date47天津大学数学系 Date48天津大学数学系 相图 Date49天津大学数学系 单摆系统的动态演示一：pendulumdyn %pendulumdyn.m function pendulumdyn g=9.8;L=1; theta0=pi/4; Y0=theta0; 0 ;%; %TSPAN=0; 10; TSPAN=0:0.05:10; T,Y=ode45(odefunp,TSPAN,Y0); Date50天津大学数学系 单摆系统的动态演示一：pendulumdyn np=length(T); figure(1) plot(T,Y(:,1);%T(1),Y(1,1), %定义小球 ball = line(color,r,Linestyle,., erasemode, xor,Markersize,30); Date51天津大学数学系 单摆系统的动态演示一：pendulumdyn for k = 1 : np set(ball,xdata,T(k),ydata,Y(k,1); drawnow; pause(0.00005); end pause Date52天津大学数学系 单摆系统的动态演示一：pendulumdyn x0=L*sin(theta0);y0=-L*cos(theta0); figure(2) plot(-0.005;0.005,0;0,color,k,linestyle,- ,linewidth,10); axis(-0.75,0.75,-1.25,0.1);axis(off);%x0,y0, Head=line(color,r,linestyle,.,erasemode,xor,m arkersize,30); body=line(0;x0,0;y0,color,b,linestyle,- ,erasemode,xor); Date53天津大学数学系 单摆系统的动态演示一：pendulumdyn for k = 1 : np theta = real(Y(k,1); theta = Y(k,1); x=L*sin(theta); y=-L*cos(theta); set(Head,xdata,x,ydata,y); set(body,xdata,0;x,ydata,0;y); drawnow; pause(0.0005); end Date54天津大学数学系 单摆系统的动态演示一：pendulumdyn %运动方程 function y = odefunp(t,y) g = 9.8; L = 0.5; c = 0.5; alpha = 1; y=y(2);-g./L.*sin(y(1)-c.*y(2).alpha; Date55天津大学数学系 Date56天津大学数学系 Date57天津大学数学系 Date58天津大学数学系 单摆系统的动态演示二：pend.m Date59天津大学数学系 非线性方程和方程组的数值解法 【例A.7.4】解下列非线性方程组，初始向 量是 。 Date60天津大学数学系 非线性方程(组)的数值解法(c1) 解：解法程序为： x0=1;1;1;-1; opt=optimset(display,iter,LargeScale,off) ; X,FVAL,EXITFLAG,OUTPUT,J = fsolve(f704,x0,opt) 结果显示为： Date61天津大学数学系 非线性方程(组)的数值解法(c2) Optimization terminated successfully: Search direction less than tolX 零点的近似数值解 X =0.18578，0.058074 ，0.77491，-0.61786。 在近似零点的目标函数值为FVAL = 0， -8.8818e-016，2.2204e-016， -2.2204e-016 EXITFLAG = 1 表明算法成功。 BACK Date62天津大学数学系 实例：一步定常迭代算法，SOR法 Step 1. 给定初始向量x(0)Rn，控制精度 。置k = 0（k为记录迭代次数的计数器） 。 Step 2. 利用一步定常迭代法计算x(k+1), x(k+1) = M x(k) + f ， M是迭代矩阵，f是常数向量。 Step 3. 如果| x(k+1) x(k)| 2 error(Too many input arguments in plegendre.m!) end Date69天津大学数学系 附录(c5) plegendre.m if nargout2 error(Too many output arguments in plegendre.m!) end switch n case 0 orthp1=1; case 1 or