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 =3
ans =3
Z =
[  237/80,  -61/80]
[  173/40, -109/40]
[ -151/40,  103/40]
[       1,       0]
[       0,       1]
x0 =00000
x =(237*k1)/80 - (61*k2)/80(173*k1)/40 - (109*k2)/40(103*k2)/40 - (151*k1)/40k1k2

 

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

syms x y z
equ1=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.x
y=root.y
z=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) t
x0=[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.x
y=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.x
y=answer.y

%%%%%%%结果%%%%%%%%%
x =
(15*exp(-t))/4 - (8*exp(-2*t))/3 - exp(-5*t)/12
y =
(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;endi=i+1;
end
x
%调用结果
x =1.00000.0000

 

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

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

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

%3-5
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];
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)+x0equ1='x.^2+y.^2=3*x*y^2';
equ2='x.^3-x.^2=y.^2-y';
ezplot(equ1)
hold on
ezplot(equ2)%%%%%%%%%%%
f=@(x)(x.^2*sin(0.1.*x+2)-3)
a=-4;
b=0
root=[]
while abs(b-a)>1e-10temp=(b+a)./2;if f(b).*f(temp)<0a=temp;elseif f(a).*f(temp)<0b=temp;elseif f(b)==0root=[root,b];elseif f(a)==0root=[root,a];elseroot=[root,temp];endbreak;end
end
root=[root,temp]%3-12
f=@(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-2equ1='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.xroot.y%3-20f=@(x)(x(1).^2-2*x(1)+x(2))g=@(t,x)(4*x(1)*x(1)+x(2)*x(2)-4)