1、M多元函数微积分常微分方程多项式矩阵线性方程组特征值 %多元函数的极限syms x ylimit(limit(x2+y2)/(sin(x)+cos(y),0),pi)limit(limit(1-cos(x2+y2)/(x2+y2)*exp(x2*y2),0),0)ans =-pi2ans =0多元函数的偏导数及全微分： clearsyms x y z dx dy dz zx zy zxx zyyz=atan(x2*y)zx=diff(z,x)zy=diff(z,y)dz=zx*dx+zy*dyzxx=diff(zx,x)zxy=diff(zx,y)z =atan(x2*y)zx =2*x*y/
2、(1+x4*y2)zy =x2/(1+x4*y2)dz =2*x*y/(1+x4*y2)*dx+x2/(1+x4*y2)*dyzxx =2*y/(1+x4*y2)-8*x4*y3/(1+x4*y2)2zxy =2*x/(1+x4*y2)-4*x5*y2/(1+x4*y2)2 zxx=diff(z,x,2)zxx =2*y/(1+x4*y2)-8*x4*y3/(1+x4*y2)2 pretty(zxx) 4 3 y x y 2 - - 8 - 4 2 4 2 2 1 + x y (1 + x y )多元函数的极值：x=fminsearch(ff,x0)x=fminunc(ff,x0) fun=i
3、nline(x(1)4+x(2)4-4*x(1)*x(2)-5);x,g=fminsearch(fun,0,0)x = 1.0000 1.0000g = -7.0000 y,h=fminunc(fun,0.1,0.1) %可以得到同样的结论多元函数积分：符号法： syms x ys=vpa(int(int(xy,x,0,1),y,1,2)s =.40546510810816438197801311546432数值法： zz=inline(x.y,x,y);s=dblquad(zz,0,1,1,2)s =0.4055内积分限为函数：先画图：x=.001:.01:3;y1=1./(2*x);y2=
4、sqrt(2*x);plot(x,y1,o,x,y2,*,2.5,-.5:.01:3)axis(-.5 3 -.5 3)legend(y1=1./(2*x),y2=sqrt(2*x),x=2.5) 确定积分限： syms x yy1=(2*x*y=1);y2=(y-sqrt(2*x)=0);x,y=solve(y1,y2,x,y)x =1/2y =1得到解： syms x yf=exp(-x2-y2);y1=1/(2*x);y2=sqrt(2*x);jfy=int(f,y,y1,y2);jfx=int(jfy,x,.5,2.5);jf2=vpa(jfx)jf2 =.12412798808725
5、833867150108282287 常微分方程求解：dsolve(eqn,var) %eqn是常微分方程，var是变量，默认是tdsolve(eqn1,eqn2,.eqnm,var)通解： s=dsolve(Dx=y,Dy=-x),y=s.y,x=s.xs = x: 1x1 sym y: 1x1 symy =C1*sin(t)+C2*cos(t)x =-C1*cos(t)+C2*sin(t)特解：dsolve(eqn,condition1,.condition n,var) y=dsolve(D2y=cos(2*x),y(0)=1,Dy(0)=0,x)simplify(y)y =-1/4*c
6、os(2*x)+5/4ans =-1/4*cos(2*x)+5/4解析解： dsolve(Dy=y-2*x/y,y(0)=1,x)ans =(2*x+1)(1/2)数值解： clearx=0:.1:4;y=sqrt(1+2*x);odefun=inline(s-2*t/s,t,s);t,s=ode45(odefun,0,4,1);plot(x,y,o-,t,s,*-)高阶微分方程式必须等价变化为一阶微分方程组：D2x-1000(1-x2)*Dx+x=0,x(0)=0,Dx(0)=1,(默认自变量为t)令：y1=x,y2=Dy1，则原微分方程组可化为： Dy1=y2,Dy2= 1000(1-y1
7、2)*y2-y1,y1(0)=0,y2(0)=1.建立程序：function dy=weifen1(t,y)dy=zeros(2,1);dy(1)=y(2);dy(2)=1000*(1-y(1)2)*y(2)-y(1);% t取值：0 3000,初值：0 1T,Y=ode15s(weifen1,0 3000,0 1);plot(T,Y(:,1),go-)求刚性微分方程： Dy1=-.01*y1-99.99*y2, Dy2=-100*y2, y1(0)=0,y2(0)=1.function f=weifen2(t,y)f=-.01*y(1)-99.99*y(2),-100*y(2);% t取值：
8、0 400,初值：2,1 t,y=ode15s(weifen2,0 400,2 1);plot(t,y,m*-)多项式的表达式与根的求解：poly2sym(p) %由多项式系数转为多项式函数polyval(p,a) %求p(a)roots(p) %求所有复数根poly(r) %由根向量r求多项式 p=1 0 -2 3;px=poly2sym(p)px =x3-2*x+3 y=polyval(p,2)y = 7 p=1 0 -2 1;q=1 0 2 3;x1=roots(p)x2=roots(q)x1 = -1.6180 1.0000 0.6180x2 = 0.5000 + 1.6583i 0.
9、5000 - 1.6583i -1.0000 多项式的四则运算：syms xconv(p1,p2) %求p1(x)与p2(x)的乘积q,r=deconv(p1,p2) %q(x)=p1(x)/p2(x)+r(x)求两多项式的和： clearp1=1 0 -2 1;p2=-1 0 1 0 -2 3;m=length(p1);n=length(p2);t=max(m,n);p1=zeros(1,t-m),p1;p2=zeros(1,t-n),p2;p=p1+p2p = -1 0 2 0 -4 4 %求两多项式的积p1=1 0 -2 1;p2=-1 0 1 0 -2 3;p3=conv(p1,p2)
10、p=poly2sym(p3)p3 = -1 0 3 -1 -4 4 4 -8 3p =-x8+3*x6-x5-4*x4+4*x3+4*x2-8*x+3求多项式的除： clearp=-1 0 3 -1 -4 4 4 -8 3;p1=1 0 -2 1;p2=1 0 -2 3;q1,r1=deconv(p,p1)q2,r2=deconv(p,p2)poly2sym(q1)poly2sym(r1)poly2sym(q2)poly2sym(r2)q1 = -1 0 1 0 -2 3r1 = 0 0 0 0 0 0 0 0 0q2 = -1 0 1 2 -2 5r2 = 0 0 0 0 0 0 -6 8
11、-12ans =-x5+x3-2*x+3ans =0ans =-x5+x3+2*x2-2*x+5ans =-6*x2+8*x-12多项式的合并与分解：syms xcollect(f) %合并同类项expand(f) %展开horner(f) %嵌套分解factor(f) %因式分解 syms x tf1=(x-1)*(x-2)*(x-3);f2=(1+x)*t+t*x;p1=collect(f1)p2=collect(f2)p1 =-6+x3-6*x2+11*xp2 =2*t*x+t clearsyms xf1=-6+x3-6*x2+11*x;f2=x5-2*x2+1;p1=horner(f1
12、)p2=horner(f2)p1 =-6+(11+(-6+x)*x)*xp2 =1+(-2+x3)*x2 clearsyms xf1=-6+x3-6*x2+11*x;f2=x5-2*x2+1;p1=factor(f1)p2=factor(f2)p1 =(x-1)*(x-2)*(x-3)p2 =(x-1)*(x4+x3+x2-x-1)有理分式的分解与合并：syms xa,b,c=residue(p,q) %将p(x)/q(x)分解为最简分式之和p,q=residue(a,b,r) %将简单分式之和合并为有理分式 clearp=1 0 0 0 0 0;q=1 -1 -2;a,b,r=residue
13、(p,q) %a表示分子a = 10.6667 0.3333b = 2 -1r = 1 1 3 5即：x5/(x2-x-2) = 10.6667/(x-2)+0.3333/(x+1)+x2+x3+3*x+5 format rataa = 32/3 1/3 将有理式分解为最简分式之和： clearp=1 0 1;q=1 -6 11 -6;a,b,r=residue(p,q)a = 5.0000 -5.0000 1.0000b = 3.0000 2.0000 1.0000r = 将有理式 分解为最简式之和： p=1;q=1 0 0 0 -1;a,b,r=residue(p,q)a = 0.2500
14、 0.0000 + 0.2500i 0.0000 - 0.2500i -0.2500 b = 1.0000 -0.0000 + 1.0000i -0.0000 - 1.0000i -1.0000 r = 将简单分式合并： cleara=32/3 1/3;b=2,-1;r=1 1 3 5;p,q=residue(a,b,r)p = 1 0 0 0 0 0q = 1 -1 -2计算行列式： D=2 -3 -1 2;1 -5 3 -4;0 2 1 -1;-5 1 3 -3;det(D)ans = -75 syms a bA=a2 a*b b2;2*a a+b 2*b;1 1 1;det(A)ans
15、=a3-3*a2*b+3*a*b2-b3克莱姆法则求解线性方程组： D=1 1 1;2 1 -1;1 -1 -1; %系数矩阵D1=2 1 1;-1 1 -1;0 -1 -1;D2=1 2 1;2 -1 -1;1 0 -1;D3=1 1 2;2 1 -1;1 -1 0;x1=det(D1)/det(D);x2=det(D2)/det(D);x3=det(D3)/det(D);x1,x2,x3x1 = 1x2 = -1x3 = 2 矩阵的生成：zeros(m,n)ones(m,n)rand(m,n) %0到1之间的均匀分布randn(m,n) %服从标准正态分布magic(n) diag(M)
16、%取上三角矩阵triu(M)tril(M)length(M)size(M)eye(n)hilb(n) %病态矩阵 a(i,j)=1/(i+j-1)pascal(n) A=magic(3)A = 8 1 6 3 5 7 4 9 2 A1=diag(A),A2=diag(A1)A1 = 8 5 2A2 = 8 0 0 0 5 0 0 0 2 triu(A),triu(A,1),triu(A,-1),ans = 8 1 6 0 5 7 0 0 2ans = 0 1 6 0 0 7 0 0 0ans = 8 1 6 3 5 7 0 9 2 A=hilb(4),A1=pascal(4)A = 1.000
17、0 0.5000 0.3333 0.2500 0.5000 0.3333 0.2500 0.2000 0.3333 0.2500 0.2000 0.1667 0.2500 0.2000 0.1667 0.1429A1 = 1 1 1 1 1 2 3 4 1 3 6 10 1 4 10 20矩阵的取块和变换：A(i,:)A(:,j)A(:) %排成一列A(i,j)A(i:j,:)A(:,i:j)A(i:j,k:l)B=reshape(A,m,n)B=rot90(A) %逆时针旋转90度B=fliplr(A)B=flipud(A)A(a,:)= % 空矩阵的一部分，会改变维数A(Aa)=b %把大
18、于a之处重新赋值B=sort(A) %把每一列的元素从小到大排列 B=A1(1:2,1:2),B1=B(:)B = 1 1 1 2B1 = 1 1 1 2 A=hilb(3),B=sort(A)A = 1.0000 0.5000 0.3333 0.5000 0.3333 0.2500 0.3333 0.2500 0.2000B = 0.3333 0.2500 0.2000 0.5000 0.3333 0.25001.0000 0.5000 0.3333矩阵的基本运算：A+(-)B %同维A+k %所有元素都加kA*B %矩阵的乘法A*k,k*A %所有元素都乘以kA.*B %同维矩阵对应元素相
19、乘A./BA/BABA. %A的转置A %A的共轭转置，在实数域内就是转置A(-1)或inv(A) %矩阵的逆AkA.kA(1/2)或sqrtm(A)sqrt(A) %矩阵中的元素开方 %求解矩阵方程：X*A=B,A*Y=BA=1 2;3 4;B=1 2;-1 0;X=B*inv(A)X1=B/AY=inv(A)*BY1=ABX = 1.0000 0 2.0000 -1.0000X1 = 1 0 2 -1Y = -3.0000 -4.0000 2.0000 3.0000Y1 = -3.0000 -4.00002.0000 3.0000 %求矩阵的范数A=1 2;3 4;A1=norm(A,1)
20、A2=norm(A,2)A3=norm(A,inf)A1 = 6A2 = 5.4650A3 = 7 %求矩阵的条件数A=1 2;3 4;A1=cond(A)A1 = 14.9330%求解线性方程组rank(A) %秩rref(A) %行的最简形null(A) %系数矩阵为A的基础解系null(A,r) %有理数形式的基础解系 A=magic(4)zhi=rank(A)zjh=rref(A)jcjx=null(A)A = 16 2 3 13 5 11 10 8 9 7 6 12 4 14 15 1zhi = 3zjh = 1 0 0 1 0 1 0 3 0 0 1 -3 0 0 0 0jcjx
21、= -0.2236 -0.6708 0.6708 0.2236 A=magic(4);B=1 2 3 4;Az=A,B,a1=rank(A),a2=rank(Az)Az = 16 2 3 13 1 5 11 10 8 2 9 7 6 12 3 4 14 15 1 4a1 = 3a2 = 4当系数矩阵的秩不等于增光矩阵的秩时，原方程组无解，但可以求出最小二乘解（X=A/b）（即A*X-b取最小值时的解）当系数矩阵的秩等于增光矩阵的秩，且秩等于方程组中方程的个数时，原方程组有唯一解，用X=A/b可求出当系数矩阵的秩等于增光矩阵的秩，且秩小于方程组中方程的个数时，原方程组有无穷多解，则X=A/b求得
22、的是一个特解jcjx=null(A) 可求出基础解系，则原方程的通解为：xt=x+k1*jcjx,k1为任意实数 X=ABWarning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.625398e-018.X = 1.0e+015 * 1.5897 4.7691 -4.7691 -1.5897一般情况下：设解得X=a,b,则可用 对原数据进行拟合特征值和特征多项式：trace(A) %迹poly(A) %矩阵A特征多项式的系数a,b=eig(A) %A的特征列向量a，对角矩
23、阵b为特征值B=orth(A) %正交化空间 A=1 2;3 4;a=trace(A),b=poly(A),c=roots(b)a = 5b = 1.0000 -5.0000 -2.0000c = 5.3723 -0.3723 AA = 1 2 3 4 c,d=eig(A)c = -0.8246 -0.4160 0.5658 -0.9094d = -0.3723 0 0 5.3723 orth(A)ans = -0.4046 -0.9145 -0.9145 0.4046 %求一个正交矩阵p，使得inv(P)*A*p=B，B为对角矩阵A=0 -1 1;-1 0 1;1 1 0;a,b=eig(A
24、)p=orth(a)B=inv(p)*A*p %验证p*inv(p) %验证a = -0.5774 -0.5216 0.6282 -0.5774 0.8048 0.1376 0.5774 0.2832 0.7658b = -2.0000 0 0 0 1.0000 0 0 0 1.0000p = -0.5774 0.7951 0.1856 -0.5774 -0.2368 -0.7814 0.5774 0.5583 -0.5958B = -2.0000 0.0000 -0.0000 0.0000 1.0000 0.0000 -0.0000 -0.0000 1.0000ans = 1.0000 -0
25、.0000 -0.0000 -0.0000 1.0000 0.00000.0000 0.0000 1.0000若不要求正交变换，则只要a可逆，则a即为所求： inv(a)ans = -0.5774 -0.5774 0.5774 -0.5216 0.8048 0.2832 0.6282 0.1376 0.7658 inv(a)*A*aans = -2.0000 -0.0000 0.0000 0.0000 1.0000 0.00000.0000 -0.0000 1.0000求一个正交变换 x=P*y,把二次型化为标准形 %二次型矩阵为：A=0 1 1 -1;1 0 -1 1;1 -1 0 1;-1
26、 1 1 0a,b=eig(A)P=orth(a)B=P*A*PP*PA = 0 1 1 -1 1 0 -1 1 1 -1 0 1 -1 1 1 0a = -0.5000 -0.0788 0.2887 -0.8127 0.5000 0.6644 -0.2887 -0.4746 0.5000 -0.7432 -0.2887 -0.3381 -0.5000 0 -0.8660 0b = -3.0000 0 0 0 0 1.0000 0 0 0 0 1.0000 0 0 0 0 1.0000P = -0.5000 0.2887 -0.8127 -0.0788 0.5000 -0.2887 -0.47
27、46 0.6644 0.5000 -0.2887 -0.3381 -0.7432 -0.5000 -0.8660 -0.0000 -0.0000B = -3.0000 0.0000 -0.0000 -0.0000 0.0000 1.0000 0.0000 0.0000 -0.0000 0.0000 1.0000 0.0000 -0.0000 0.0000 0.0000 1.0000ans = 1.0000 -0.0000 0.0000 0.0000 -0.0000 1.0000 0.0000 0.0000 0.0000 0.0000 1.0000 0.00000.0000 0.0000 0.0000 1.0000即得标准形为：将P的列向量移动可相应的交换标准型的系数
