Matlab实验（二）：典型操作与代码示例

Matlab 实验 （二）

1.求积分

f=@(x,y)(x.*cos(x+y.*y));dblquad(f,0,pi,0,2*pi)   ans =   -3.4267

2.求常微分方程的 数值解

syms y;dsolve('D2y+2*Dy+y=0','y(0)=0','Dy(0)=1')  ans =t*exp(-t)

3.试求下面齐次方程的基础解系。

%调用代码A=[+6  +1 +4 -7 -3;    -2 -7 -8 +6 0;    -4 +5 +1 -6 +8;    -34 +36 +9 -21 49;     -26 -12 -27 27 17;];B=zeros(5,1)；rank(A)rank([A,B])%判断两者秩Z=null(sym(A))%Ax=0，此时已经是解了x0=sym(pinv(A))*B%可以不用syms k1 k2;%求出结果x=k1*Z(:,1)+k2*Z(:,2)+x0%验证A*x-B

%调用结果ans =     3ans =     3Z =[  237/80,  -61/80][  173/40, -109/40][ -151/40,  103/40][       1,       0][       0,       1]x0 = 0 0 0 0 0x =   (237*k1)/80 - (61*k2)/80 (173*k1)/40 - (109*k2)/40 (103*k2)/40 - (151*k1)/40                        k1                        k2

4.试求下面的代数方程的全部的根

syms x y zequ1=x^2*y^2-x*y*z-4*x^2*y*z^2==x*z^2;equ2=x*y^3-2*y*z^2==3*x^3*z^2+4*x*z*y^2;equ3=y*x^2-7*x*y+3*x*z^2==x^4*z*y;equs=[equ1,equ2,equ3]root=solve(equs);x=root.xy=root.yz=root.z%由于结果是符号解，很长，这里就不展现。

5.Lorenz方程是研究混沌问题的著名非线性微分方程，其数学形式为

其中， ，且其初值为 。试求出其数值解，绘制三维空间曲线，并绘制Lorenz方程解在两平面上的投影。

function dx=func_2_5(t,x)beta=8/3;theta=10;gamma=28;dx=[-beta*x(1)+x(2)*x(3);    -theta*x(2)+theta*x(3);    -x(1)*x(2)+gamma*x(2)-x(3)];end

syms x(t) tx0=[0 0 10^(-3)];[t,x]=ode45('func_2_5',[0,200],x0); figure(1);plot(t,x)grid on;axis auto;xlabel('t');ylabel('x');title('二维Lorenz方程'); figure(2);plot3(x(:,1),x(:,2),x(:,3))grid on;axis auto;xlabel('x(1)');ylabel('x(2)');zlabel('x(3)');title('三维Lorenz方程'); figure(3);comet3(x(:,1),x(:,2),x(:,3))grid on;axis auto;xlabel('x(1)')ylabel('x(2)');zlabel('x(3)');title('三维Lorenz方程'); view(0,90)view(0,0)

6.试求出下面微分方程组的解析解，并和数值解比较。

%解析解equ1='D2x=-2*x-3*Dx-exp(-5*t)';equ2='D2y=2*x-3*y-4*Dx-4*Dy-sin(t)';equ3='x(0)=1';equ4='Dx(0)=2';equ5='y(0)=3';equ6='Dy(0)=4';answer=dsolve(equ1,equ2,equ3,equ4,equ5,equ6);x=answer.xy=answer.y  %第二种调用方法syms y(t) x(t)eqns = [diff(x,t,2) == -2*x-3*diff(x,t)-exp(-5*t), diff(y,t,2) == 2*x-3*y-4*diff(x,t)-4*diff(y,t)-sin(t)];Dx=diff(x,t);Dy=diff(y,t);conds = [x(0)==1,Dx(0)==2,y(0)==3,Dy(0)==4];answer=dsolve(eqns,conds);x=answer.xy=answer.y

%%%%%%%结果%%%%%%%%%x =(15*exp(-t))/4 - (8*exp(-2*t))/3 - exp(-5*t)/12y =(80*exp(-2*t))/3 - (207*exp(-t))/16 - (107*exp(-3*t))/10 - (11*exp(-5*t))/48 + (45*t*exp(-t))/4 + (5^(1/2)*cos(t + atan(1/2)))/10

%%%%%%%%%%%数值解%%%%%%%%function dx=func_2_6(t,x)dx=[x(2);    -2*x(1)-3*x(2)+exp(-5*t);    x(4); 2*x(1)-3*x(3)-4*x(2)-4*x(4)-sin(t);    ];end

%调用tspan=[-5 5];xy0 =[1;2;3;4];[t x]=ode45(@func_2_6,tspan,xy0);%作图plot(x(:,2),x(:,4))

7.求解下面的 最优化 问题

function sx=func_2_7_1(x)sx=x(1)^2-2*x(1)+x(2);end

function [sx se]=func_2_7_2(x)se=[];sx=[4*x(1)^2+x(2)^2-4];end

x0=[5;5];A=[];B=[];Aeq=[];Beq=[];xm=[0;0];xM=[;];ff=optimset;ff.Tolx=1e-10;ff.TolFun=1e-20;   %x=fmincon('func_2_7_1', x0,A,B,Aeq,Beq,xm,xM,'func_2_7_2',ff)i=1;while(1)    [x,a,b]=fmincon('func_2_7_1', x0,A,B,Aeq,Beq,xm,xM,'func_2_7_2',ff)    if b>0,break;end    i=i+1;endx%调用结果x =    1.0000    0.0000

—————————————————————————————————————————————————————————————————————————————————

说明：以上习题来自《控制系统 计算 机辅助设计 —— MATLAB语言与应用( 第3版)》的第三章，薛定宇老师的著作，解答仅供参考，如有错误，请指正。

以下是第三章习题的一些解答，仅供参考。

%3-5A=[6 1 4 -7 -3;    -2 -7 -8 6 0;    -4 5 1 -6 8;    -34 36 9 -21 49;    -26 -12 -27 27 17];rank(A)length=size(A)rank(A+zeros(length(2),1))syms k1 k2;x0=pinv(sym(A))*zeros(length(2),1)z=null(sym(A))x=k1*z(:,1)+k2*z(:,2)+x0 equ1='x.^2+y.^2=3*x*y^2';equ2='x.^3-x.^2=y.^2-y';ezplot(equ1)hold onezplot(equ2) %%%%%%%%%%%f=@(x)(x.^2*sin(0.1.*x+2)-3)a=-4;b=0root=[]while abs(b-a)>1e-10    temp=(b+a)./2;    if f(b).*f(temp)<0        a=temp;    elseif f(a).*f(temp)<0        b=temp;    else        if f(b)==0            root=[root,b];        elseif f(a)==0            root=[root,a];        else            root=[root,temp];        end        break;    endendroot=[root,temp]            %3-12f=@(t,x)([-8/3*x(1)+x(2)*x(3);-10*x(2)+10*x(3);-x(1)*x(2)+28*x(2)-x(3)])[t,y]=ode45(f,[-100 100],[0 0 1e-3]);plot(t,y) comet3(y(:,1),y(:,2),y(:,3))  %3-15-2 equ1='D2x+D2y+x+y=0'; equ2='2*D2x-D2y-x+y=sin(t)'; root=dsolve(equ1,equ2,'x(0)=2','Dx(0)=-1','y(0)=1','Dy(0)=-1') root.x root.y  %3-20 f=@(x)(x(1).^2-2*x(1)+x(2)) g=@(t,x)(4*x(1)*x(1)+x(2)*x(2)-4)

免责声明：本文系网络转载或改编，未找到原创作者，版权归原作者所有。如涉及版权，请联系删