MATLAB的科学计算.ppt
MATLAB的科学计算,5.1解析解与数值解5.2数值线性代数问题求解特殊矩阵输入zeros(m,n)ones(m,n)eye(m,n)rand(m,n)5.2.2 矩阵的特征参数运算1、行列式det(a)a=1 2 3;4 5 6;7 8 0;det(a)2、迹trace(a)对角线元素和3、秩rank(a)线性无关4、范数norm(a,选项）norm(a),norm(a,2),norm(a,1),norm(a,inf),norm(a,fro),5.2数值线性代数问题求解,5、特征多项式、特征方程、特征根构造矩阵sI-A求出其行列式,可得到多项式：C(s)=det(sI-A)=sn+c1sn-1+.+cn-1s+cnC(s)为A的特征多项式，ci为多项式系数c=poly(a)C(s)=s3-6s2-72s-27令 C(s)0构成的方程为特征方程，其根为矩阵的特征根。eig(a)多项式方程求根函数roots(c),5.2数值线性代数问题求解,6、多项式及多项式矩阵的求值polyval(aa,x)polyvalm(aa,A)aa:降幂排列多项式向量。x:给定变量(标量），数组规则运算。A:给定矩阵，矩阵规则运算。a=1 2 3;4 5 6;7 8 0 aa=poly(a)b=polyval(aa,2)c=polyvalm(aa,a),5.2数值线性代数问题求解,多项式处理函数,5.2数值线性代数问题求解,矩阵的相似变换与分解1、三角形分解LU分解任何一个方阵可以表示成两个三角矩阵的乘积，其中一个是换位的下三角阵，一个是上三角阵。lu函数可得到分解后的两个三角阵。a=1 2 3;4 5 6;7 8 0;l,u=lu(a)l*u,det(a),det(l)*det(u),5.2数值线性代数问题求解,2、正交分解QR分解将任何矩阵表示成一个正交矩阵和一个上三角矩阵的乘积。q,r=qr()a=1 2 3;4 5 6;7 8 9;10 11 12q,r=qr(a)q*r3、奇异值分解u,s,v=svd(a)4、特征值特征向量x,d=eig(a)d:特征值，x:每一列是一个特征向量。,5.3数值微积分问题,Matlab中有一类函数不对数值矩阵工作，而对数学函数工作。数值积分、非线性方程求解、优化、微分方程求解。5.3.1数值差分运算dy=diff(y)y由yi,i=1 2.n组成,diff处理后得出一个新向量yi+1-yi,i=1 2.nv=vander(1:6)%万达摩方阵diff(v)对每一列进行差分运算。列数不变，行数减1 dx,dy=gradient(a)二维差分运算,5.3数值微积分问题,5.3.2数值积分函数定积分的数值方法很多，其基本思想是将整个积分空间分割为若干子空间，子空间可积，整体可积。Matlab基于这样的思想采用自适应变步长方法给出quad()函数求定积分。,5.3数值微积分问题,编写humps.m文件function y=humps(x)y=1./(x-.3).2+0.01)+1./(x-.9).2+0.04)-6;绘图x=-1:0.01:2;plot(x,humps(x)01积分q=quad(humps,0,1)一般调用格式：y,n=quad(F,a,b,tol)y:积分值、n：被积函数调用次数、F：被积函数，一般用一个.m文件表示。a、b 上下限、tol:变步长积分误差限，默认1e-3。,5.3数值微积分问题,常用一元函数积分指令p93:表53quidquid8quid1trapzsumfnintquid与quid8误差线1e-6,新引入quid1代替。除直接用.m文件描述函数外，还可使用inline()函数定义一个函数例：f=inline(1/sqrt(2*pi)*exp(-x.2/2),x),y,kk=quad(f,-8,8)避免建立不必要的文件双重积分：dblquad(函数名，xm,Xm,ym,Ym,yol)I=,5.4常微分方程的数值解法,一般常微分方程的数值解法：只包含一个自变量的微分方程（ODE问题）：常微分方程。初值问题、边值问题假定一阶常微分方程xi=fi(t,x),i=1,2,3.,nx 状态变量 xi 构成的向量，x=x1,x2,.xn T状态向量。n为系统阶次fi(.)为任意非线性函数t时间变量设定初值x(o),求解常微分方程组。表5-4 p95t,x=ode23(方程函数名，t0,tf,x0,选项)t,x=ode45(方程函数名，t0,tf,x0,选项)选项获取设置：odeget()odeset()方程函数名:描述系统状态方程的M函数名称，用引号括起来。其编写格式固定function xdot=方程函数名(t,x)t0 tf：用户指定起始和终止计算时间。x0：系统初始状态变量值。返回t:求解时间变量；x:各时刻状态变量列向量构成的矩阵转置。,5.4常微分方程的数值解法,x1(t)=-8x1(t)/3+x2(t)x3(t)x2(t)=-10 x2(t)+10 x3(t)x3(t)=-x1(t)x2(t)+28x2(t)-x3(t)x0=0;0;1e-10动态函数模型lorenzeq.m用ode45对lorenzeq描述系统状态方程进行数值求解。function xdot=lorenzeq(t,x)xdot=-8/3*x(1)+x(2)*x(3);-10*x(2)+10*x(3);-x(1)*x(2)+28*x(2)-x(3);t_final=100;x0=0;0;1e-10;t,x=ode45(lorenzeq,0,t_final,x0);figure(1);set(gcf,position,7,225,390,270)%获取当前窗口句柄并设置plot(t,x)figure(2);set(gcf,position,405,225,390,270)plot3(x(:,1),x(:,2),x(:,3);axis(10 40-20 20-20 20),5.4常微分方程的数值解法,常微分方程组的变换与技巧ode函数只能解决一阶常微分方各组问题，对于高阶常微分问题，应首先将其变换为一阶常微分方程问题。1.单个高阶常微分方程处理方法一个高阶常微分方程的一般形式y(n)=f(t,y,y,.,y(n-1)其已知输出变量y(t)的各阶导数初始值为y(o),y(0),.y(n-1)(0),则可以选择一组状态变量x1=y,x2=y.xn=y(n-1),将原高阶方程变换成如下的一阶方程组。x1=x2x2=x3.xn=f(t,x1,x2,.xn)初值：x1(0)=y(0),x2(0)=y(0),.xn(0)=y(n-1)(0).,5.4常微分方程的数值解法,2.高阶常微分方程组的变换x(m)=f(t,x,x.x(m-1),y,y,.y(n-1)y(n)=g(t,x,x.x(m-1),y,y,.y(n-1)令：x1=x,x2=x,.,xm=x(m-1),x(m+1)=y,x(m+2)=y,.,x(m+n)=y(n-1)x1=x2x2=x3.xm=f(t,x1,x2,.,xm+n)xm+1=xm+2xm+2=xm+3.xm+n=g(t,x1,x2,.,xm+n),5.4常微分方程的数值解法,例：二阶微分方程x+(x2-1)x+x=0function xdot=difeq(t,x)xdot=zeros(2,1);xdot(1)=x(1).*(1-x(2).2)-x(2);xdot(2)=x(1);t0=0;tf=20;x0=0;0.25;t,x=ode23(difeq,t0,tf,x0);plot(t,x),5.4常微分方程的数值解法,5-3Mx=Mg-kxx=g-k/mxfunction dx=sink(t,x)m=1;g=9.81;k=10;dx=x(2);g-k/m*x(2);x0=0;1;t,x=ode45(sink,0,1,x0);subplot(2,1,1),plot(t,x(:,1);subplot(2,1,2),plot(t,x(:,2);,5.5 数值插值与统计分析,一维数据的插值拟合f(x)为一维函数，由已知点求出其它点的函数值。1、一维插值函数interp1()y1=interp1(x,y,x1,方法）x,y 已知点数据；x1一组新插值点；y1函数值；方法：linear线性；cubic三次；spline样条型2、多项式拟合函数polyfit()p=polyfit(x,y,n)x y 已知一组自变量和函数值数据；n:预期多项式的阶次，p：返回多项式系数,5.5 数值插值与统计分析,p=0,1.1,2.1,2.8,4.2,5,6.1,6.9,8.1,9,9.9;u=10,11,13,14,17,18,22,24,29,34,39;a=polyfit(p,u,3);disp(a)p1=0:0.01:10;u1=polyval(a,p1);plot(p1,u1,p,u,o)gtext(datestr(today)一维插值点比较x1=0.6,3.2,7.5,8.5;y1=interp1(p,u,x1);plot(p1,u1,p,u,o,x1,y1,p)二维数据插值拟合interp2(),5.5 数值插值与统计分析,数据分析与统计处理表55不论对行向量还是列向量，函数均按向量长度处理，对矩阵则按列处理。x,i=max(v),各列最大值置于经x,其位置行号置于i.a=magic(8)x,i=max(a)max(max(a),
