东北大学MATLAB实验答案

东北大学MATLAB实验答案 东北大学《MATLAB语言与应用》实验 东北大学MATLAB实验课习答案 第一部分MATLAB语言编程、科学绘图与基本问题求解 2、 >> A=[1,2,3,4;4,3,2,1;2,3,4,1;3,2,4,1] A = 1 2 3 4 4 3 2 1 2 3 4 1 3 2 4 1 B=[1+4j,2+3j,3+2j,4+1j;4+1j,3+2j,2+3j,1+4j;2+3j,3+2j,4+1j,1+4j;3+2j,2+3j,4+1j,1+4j] B = 1.0000 + 4.0000i 2.0000 + 3.0000i 3.0000 + 2.0000i 4.0000 + 1.0000i 4.0000 + 1.0000i 3.0000 + 2.0000i 2.0000 + 3.0000i 1.0000 + 4.0000i 2.0000 + 3.0000i 3.0000 + 2.0000i 4.0000 + 1.0000i 1.0000 + 4.0000i 3.0000 + 2.0000i 2.0000 + 3.0000i 4.0000 + 1.0000i 1.0000 + 4.0000i >> A(5,6)=5 A = 1 2 3 4 0 0 4 3 2 1 0 0 2 3 4 1 0 0 3 2 4 1 0 0 0 0 0 0 0 5 3、 A=magic(8) A = 64 2 3 61 60 6 7 57 9 55 54 12 13 51 50 16 17 47 46 20 21 43 42 24 40 26 27 37 36 30 31 33 32 34 35 29 28 38 39 25 41 23 22 44 45 19 18 48 49 15 14 52 53 11 10 56 8 58 59 5 4 62 63 1 >> B=A(2:2:end,:) B = 9 55 54 12 13 51 50 16 40 26 27 37 36 30 31 33 41 23 22 44 45 19 18 48 8 58 59 5 4 62 63 1 4( i=0:63;s=sum(2.^i) s = 1.8447e+019 5、 (1) >> z=sin(1./t); Warning: Divide by zero. >> plot(t,z) 1 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1-1-0.8-0.6-0.4-0.200.20.40.60.81 (2) >> t=[-pi:0.05:-1.8,-1.799:.001:-1.2,-1.2:0.05:1.2,1.201:0.001:1.8,1.81:0.05:pi]; >> y=sin(tan(t))-tan(sin(t)); >> plot(t,y) 3 2 1 0 -1 -2 -3-4-3-2-101234 6. >> xx=[-2:.1:-1.2,-1,1:0.02:-0.9,-0.8:0.1:0.8,0.9:0.02:1.1,1.2:0.1:2]; >> yy=[-1:0.1:-0.2,-0.1:0.02:0.1,0.2:.1:1]; >> [x,y]=meshgrid(xx,yy); >> z=1./(sqrt((1-x).^2+y.^2))+1./(sqrt((1+x).^2+y.^2)); >>surf(x,y,z),shading flat;zlim([0,15]) 15 10 5 0 1 0.52 100-0.5-1 -1-2 三视图 >> surf(x,y,z),shading flat;zlim([0,15]) >> xx=[-2:.1:-1.2,-1,1:0.02:-0.9,-0.8:0.1:0.8,0.9:0.02:1.1,1.2:0.1:2]; >> yy=[-1:0.1:-0.2,-0.1:0.02:0.1,0.2:.1:1]; >> [x,y]=meshgrid(xx,yy); >> subplot(224),surf(x,y,z) >> subplot(221),surf(x,y,z),view(0,90); >> subplot(222),surf(x,y,z),view(90,0); >> subplot(223),surf(x,y,z),view(0,0); 160 0.540 0 20-0.5 -10-2-1012-1-0.500.51 60 100 40 50 20012000-1-2-2-1012 7. (1) >> syms x;f=(3.^x+9.^x)^(1./x);L=limit(f,x,inf) L = 9 (2) >>syms x y;f=x*y/(sqrt(x*y+1)-1);L1=limit(limit(f,x,0),y,0) L1 = 2 (3) >> syms x y; >> f=(1-cos(x^2+y^2))/((x^2+y^2)*exp(x^2+y^2)); >> L=limit(limit(f,x,0),y,0) L = 0 8( 先建立M文件: function result=paradiff(y,x,t,n) if mod(n,1)~=0,error('n should positive integer,please correct') else if n==1,result=diff(y,t)/diff(x,t); else,result=diff(paradiff(y,x,t,n-1),t)/diff(x,t); end,end 然后调用函数: >> syms t;x=log(cos(t));y=cos(t)-t*sin(t); >> f=paradiff(y,x,t,1); >> [n,d]=numden(f); >> F=simple(n)/simple(d) F = (2*sin(t)+t*cos(t))*cos(t)/sin(t) >> syms t;x=log(cos(t));y=cos(t)-t*sin(t); >> f=paradiff(y,x,t,1); >> syms t;x=log(cos(t));y=cos(t)-t*sin(t); >> f=paradiff(y,x,t,2); >> [n,d]=numden(f); >> F=simple(n)/simple(d) F = -cos(t)*(3*cos(t)^2*sin(t)+cos(t)^3*t-2*sin(t)-2*t*cos(t))/sin(t)^3 >> subs(F,t,pi/3) ans = 1.5387 9. >> syms x y t; >> f=exp(-t^2); >> I=simple(int(f,t,0,x*y)) I = 1/2*pi^(1/2)*erf(x*y) >> F=x/y*diff(I,x,2)-2*diff(diff(I,x),y)+diff(I,y,2) F = 2*x^2*y^2*exp(-x^2*y^2)-2*exp(-x^2*y^2)-2*x^3*y*exp(-x^2*y^2) >> K=simple(F) K = -2*exp(-x^2*y^2)*(-x^2*y^2+1+x^3*y) 10. (1) >> syms n; >> S=symsum(1/((2*n)^2-1),n,1,inf) S = 1/2 (2) >> syms k n >> limit(n*symsum(1/(n^2+k*pi),k,1,n),n,inf) ans = 1 11. (1)>> syms t; >> syms a positive; >> x=a*(cos(t)+t*sin(t)); >> y=a*(sin(t)-t*cos(t)); >> I=int((x^2+y^2)*sqrt(diff(x,t)^2+diff(y,t)^2),t,0,2*pi) I = 2*a^3*pi^2+4*a^3*pi^4 (2) >> syms t; >> syms a b c positive; >> x=c/a*cos(t); >> y=c/b*sin(t); >> F=[y*x^3+exp(y),x*y^3+x*exp(y)-2*y]; >> ds=[diff(x,t);diff(y,t)]; >> I=int(F*ds,t,pi,0) I = 2/15*c*(-2*c^4+15*b^4)/a/b^4 12. 首先编写M程序: function A=vander(v) n=length(v);v=v(:);A=sym(ones(n)); for j=n-1:-1:1,A(:,j)=v.*A(:,j+1);end >> syms a b c d e; >> A=[a,b,c,d,e]; >> V=vander(A) V = [ a^4, a^3, a^2, a, 1] [ b^4, b^3, b^2, b, 1] [ c^4, c^3, c^2, c, 1] [ d^4, d^3, d^2, d, 1] [ e^4, e^3, e^2, e, 1] >> det(V),simple(ans) ans = (c-d)*(b-d)*(b-c)*(a-d)*(a-c)*(a-b)*(-d+e)*(e-c)*(e-b)*(e-a) 13. >> A=[-2,0.5,-0.5,0.5;0,-1.5,0.5,-0.5;2,0.5,-4.5,0.5;2,1,-2,-2]; >> [V J]=jordan(sym(A)) V = [ 0, 1/2, 1/2, -1/4] [ 0, 0, 1/2, 1] [ 1/4, 1/2, 1/2, -1/4] [ 1/4, 1/2, 1, -1/4] J = [ -4, 0, 0, 0] [ 0, -2, 1, 0] [ 0, 0, -2, 1] [ 0, 0, 0, -2] 14. 先编写M文件: function X=lyap(A,B,C) if nargin==2,C=B;B=A';end [nr,nc]=size(C);A0=kron(A,eye(nc))+kron(eye(nr),B'); try C1=C'; x0=-inv(A0)*C1(:);X=reshape(x0,nc,nr)'; catch,error('singular matrix found.'),end 数值解为: >> A=[3,-6,-4,0,5;1,4,2,-2,4;-6,3,-6,7,3;-13,10,0,-11,0;0,4,0,3,4]; >> B=[3,-2,1;-2,-9,2;-2,-1,9]; >> C=[-2,1,-1;4,1,2;5,-6,1;6,-4,-4;-6,6,-3]; >> X=lyap(A,B,C) X = -4.0569 -14.5128 1.5653 0.0356 25.0743 -2.7408 9.4886 25.9323 -4.4177 2.6969 21.6450 -2.8851 7.7229 31.9100 -3.7634 >> norm(A*X+X*B+C) ans = 3.9870e-013 解析解为: >> X=lyap(sym(A),B,C) X = [ -434641749950/107136516451, -4664546747350/321409549353, 503105815912/321409549353] [ 3809507498/107136516451, 8059112319373/321409549353, -880921527508/321409549353] [ 1016580400173/107136516451, 8334897743767/321409549353, -1419901706449/321409549353] [ 288938859984/107136516451, 6956912657222/321409549353, -927293592476/321409549353] [ 827401644798/107136516451, 10256166034813/321409549353, -1209595497577/321409549353] >> A*X+X*B+C ans = [ 0, 0, 0] [ 0, 0, 0] [ 0, 0, 0] [ 0, 0, 0] [ 0, 0, 0] 15. (1) >> A=[-4.5,0,0.5,-1.5;-0.5,-4,0.5,-0.5;1.5,1,-2.5,1.5;0,-1,-1,-3]; >> A=sym(A);syms t; >> expm(A*t) ans = [ 1/2*exp(-5*t)+1/2*exp(-3*t)-1/2*t*exp(-3*t)+1/2*t^2*exp(-3*t), t*exp(-3*t)-1/2*exp(-3*t)+1/2*exp(-5*t), 1/2*t^2*exp(-3*t)+1/2*t*exp(-3*t), 1/2*t^2*exp(-3*t)-1/2*t*exp(-3*t)-1/2*exp(-3*t)+1/2*exp(-5*t)] [ -1/2*exp(-3*t)+1/2*exp(-5*t)+1/2*t*exp(-3*t), 1/2*exp(-5*t)+1/2*exp(-3*t), 1/2*t*exp(-3*t), -1/2*exp(-3*t)+1/2*exp(-5*t)+1/2*t*exp(-3*t)] [ 1/2*exp(-3*t)-1/2*exp(-5*t)+1/2*t*exp(-3*t), 1/2*exp(-3*t)-1/2*exp(-5*t), 1/2*t*exp(-3*t)+exp(-3*t), 1/2*exp(-3*t)-1/2*exp(-5*t)+1/2*t*exp(-3*t)] [ -1/2*t^2*exp(-3*t), -t*exp(-3*t), -t*exp(-3*t)-1/2*t^2*exp(-3*t), exp(-3*t)-1/2*t^2*exp(-3*t)] (2) 编写M程序 function F=funm(A,fun,x) [V,J]=jordan(A); v1=[0,diag(J,1)']; v2=[find(v1==0),length(v1)+1]; for i=1:length(v2)-1 v_lambda(i)=J(v2(i),v2(i));v_n(i)=v2(i+1)-v2(i); end m=length(v_lambda);F=sym([]); for i=1:m J1=J(v2(i):v2(i)+v_n(i)-1,v2(i):v2(i)+v_n(i)-1); fJ=funJ(J1,fun,x);F=diagm(F,fJ); end F=V*F*inv(V); function fJ=funJ(J,fun,x) lam=J(1,1); f1=fun; fJ=subs(fun,x,lam)*eye(size(J)); H=diag(diag(J,1),1); H1=H; for i=2:length(J) f1=diff(f1,x); a1=subs(f1,x,lam); fJ=fJ+a1*H1; H1=H1*H/i; end function A=diagm(A1,A2) A=A1; [n,m]=size(A); [n1,m1]=size(A2); A(n+1:n+n1,m+1:m+m1)=A2; >> A=[-4.5,0,0.5,-1.5;-0.5,-4,0.5,-0.5;1.5,1,-2.5,1.5;0,-1,-1,-3]; >> syms x t; >> A1=funm(sym(A),sin(x*t),x) A1 = [ -1/2*sin(5*t)+1/2*sin(3*t)*t^2-1/2*cos(3*t)*t-1/2*sin(3*t), -1/2*sin(5*t)+cos(3*t)*t+1/2*sin(3*t), 1/2*cos(3*t)*t+1/2*sin(3*t)*t^2, -1/2*sin(5*t)+1/2*sin(3*t)+1/2*sin(3*t)*t^2-1/2*cos(3*t)*t] [ -1/2*sin(5*t)+1/2*cos(3*t)*t+1/2*sin(3*t), -1/2*sin(5*t)-1/2*sin(3*t), 1/2*cos(3*t)*t, -1/2*sin(5*t)+1/2*cos(3*t)*t+1/2*sin(3*t)] [ 1/2*sin(5*t)+1/2*cos(3*t)*t-1/2*sin(3*t), 1/2*sin(5*t)-1/2*sin(3*t), -sin(3*t)+1/2*cos(3*t)*t, 1/2*sin(5*t)+1/2*cos(3*t)*t-1/2*sin(3*t)] [ -1/2*sin(3*t)*t^2, -cos(3*t)*t, -cos(3*t)*t-1/2*sin(3*t)*t^2, -sin(3*t)-1/2*sin(3*t)*t^2] (3) >> A=[-4.5,0,0.5,-1.5;-0.5,-4,0.5,-0.5;1.5,1,-2.5,1.5;0,-1,-1,-3]; >> syms x t; >> A1=funm(sym(A),exp(x*t)*sin(x^2*exp(x*t)*t),x) A1 = [ 1/2*exp(-5*t)*sin(25*exp(-5*t)*t)+1/2*t^2*exp(-3*t)*sin(9*exp(-3*t)*t)+t*exp(-3*t)*cos(9*exp(-3*t)*t)*(-6*exp(-3*t)*t+9*t^2*exp(-3*t))-1/2*exp(-3*t)*sin(9*exp(-3*t)*t)*(-6*exp(-3*t)*t+9*t^2*exp(-3*t))^2+1/2*exp(-3*t)*cos(9*exp(-3*t)*t)*(2*exp(-3*t)*t-12*t^2*exp(-3*t)+9*t^3*exp(-3*t))-1/2*t*exp(-3*t)*sin(9*exp(-3*t)*t)-1/2*exp(-3*t)*cos(9*exp(-3*t)*t)*(-6*exp(-3*t)*t+9*t^2*exp(-3*t))+1/2*exp(-3*t)*sin(9*exp(-3*t)*t), 1/2*exp(-5*t)*sin(25*exp(-5*t)*t)+t*exp(-3*t)*sin(9*exp(-3*t)*t)+exp(-3*t)*cos(9*exp(-3*t)*t)*(-6*exp(-3*t)*t+9*t^2*exp(-3*t))-1/2*exp(-3*t)*sin(9*exp(-3*t)*t), 1/2*t*exp(-3*t)*sin(9*exp(-3*t)*t)+1/2*exp(-3*t)*cos(9*exp(-3*t)*t)*(-6*exp(-3*t)*t+9*t^2*exp(-3*t))+1/2*t^2*exp(-3*t)*sin(9*exp(-3*t)*t)+t*exp(-3*t)*cos(9*exp(-3*t)*t)*(-6*exp(-3*t)*t+9*t^2*exp(-3*t))-1/2*exp(-3*t)*sin(9*exp(-3*t)*t)*(-6*exp(-3*t)*t+9*t^2*exp(-3*t))^2+1/2*exp(-3*t)*cos(9*exp(-3*t)*t)*(2*exp(-3*t)*t-12*t^2*exp(-3*t)+9*t^3*exp(-3*t)), 1/2*exp(-5*t)*sin(25*exp(-5*t)*t)-1/2*exp(-3*t)*sin(9*exp(-3*t)*t)+1/2*t^2*exp(-3*t)*sin(9*exp(-3*t)*t)+t*exp(-3*t)*cos(9*exp(-3*t)*t)*(-6*exp(-3*t)*t+9*t^2*exp(-3*t))-1/2*exp(-3*t)*sin(9*exp(-3*t)*t)*(-6*exp(-3*t)*t+9*t^2*exp(-3*t))^2+1/2*exp(-3*t)*cos(9*exp(-3*t)*t)*(2*exp(-3*t)*t-12*t^2*exp(-3*t)+9*t^3*exp(-3*t))-1/2*t*exp(-3*t)*sin(9*exp(-3*t)*t)-1/2*exp(-3*t)*cos(9*exp(-3*t)*t)*(-6*exp(-3*t)*t+9*t^2*exp(-3*t))] [ 1/2*exp(-5*t)*sin(25*exp(-5*t)*t)+1/2*t*exp(-3*t)*sin(9*exp(-3*t)*t)+1/2*exp(-3*t)*cos(9*exp(-3*t)*t)*(-6*exp(-3*t)*t+9*t^2*exp(-3*t))-1/2*exp(-3*t)*sin(9*exp(-3*t)*t), 1/2*exp(-5*t)*sin(25*exp(-5*t)*t)+1/2*exp(-3*t)*sin(9*exp(-3*t)*t), 1/2*t*exp(-3*t)*sin(9*exp(-3*t)*t)+1/2*exp(-3*t)*cos(9*exp(-3*t)*t)*(-6*exp(-3*t)*t+9*t^2*exp(-3*t)), 1/2*exp(-5*t)*sin(25*exp(-5*t)*t)+1/2*t*exp(-3*t)*sin(9*exp(-3*t)*t)+1/2*exp(-3*t) *cos(9*exp(-3*t)*t)*(-6*exp(-3*t)*t+9*t^2*exp(-3*t))-1/2*exp(-3*t)*sin(9*exp(-3*t)*t)] [ -1/2*exp(-5*t)*sin(25*exp(-5*t)*t)+1/2*t*exp(-3*t)*sin(9*exp(-3*t)*t)+1/2*exp(-3*t)*cos(9*exp(-3*t)*t)*(-6*exp(-3*t)*t+9*t^2*exp(-3*t))+1/2*exp(-3*t)*sin(9*exp(-3*t)*t), -1/2*exp(-5*t)*sin(25*exp(-5*t)*t)+1/2*exp(-3*t)*sin(9*exp(-3*t)*t), exp(-3*t)*sin(9*exp(-3*t)*t)+1/2*t*exp(-3*t)*sin(9*exp(-3*t)*t)+1/2*exp(-3*t)*cos(9*exp(-3*t)*t)*(-6*exp(-3*t)*t+9*t^2*exp(-3*t)), -1/2*exp(-5*t)*sin(25*exp(-5*t)*t)+1/2*t*exp(-3*t)*sin(9*exp(-3*t)*t)+1/2*exp(-3*t)*cos(9*exp(-3*t)*t)*(-6*exp(-3*t)*t+9*t^2*exp(-3*t))+1/2*exp(-3*t)*sin(9*exp(-3*t)*t)] [ -1/2*t^2*exp(-3*t)*sin(9*exp(-3*t)*t)-t*exp(-3*t)*cos(9*exp(-3*t)*t)*(-6*exp(-3*t)*t+9*t^2*exp(-3*t))+1/2*exp(-3*t)*sin(9*exp(-3*t)*t)*(-6*exp(-3*t)*t+9*t^2*exp(-3*t))^2-1/2*exp(-3*t)*cos(9*exp(-3*t)*t)*(2*exp(-3*t)*t-12*t^2*exp(-3*t)+9*t^3*exp(-3*t)), -t*exp(-3*t)*sin(9*exp(-3*t)*t)-exp(-3*t)*cos(9*exp(-3*t)*t)*(-6*exp(-3*t)*t+9*t^2*exp(-3*t)), -t*exp(-3*t)*sin(9*exp(-3*t)*t)-exp(-3*t)*cos(9*exp(-3*t)*t)*(-6*exp(-3*t)*t+9*t^2*exp(-3*t))-1/2*t^2*exp(-3*t)*sin(9*exp(-3*t)*t)-t*exp(-3*t)*cos(9*exp(-3*t)*t)*(-6*exp(-3*t)*t+9*t^2*exp(-3*t))+1/2*exp(-3*t)*sin(9*exp(-3*t)*t)*(-6*exp(-3*t)*t+9*t^2*exp(-3*t))^2-1/2*exp(-3*t)*cos(9*exp(-3*t)*t)*(2*exp(-3*t)*t-12*t^2*exp(-3*t)+9*t^3*exp(-3*t)), exp(-3*t)*sin(9*exp(-3*t)*t)-1/2*t^2*exp(-3*t)*sin(9*exp(-3*t)*t)-t*exp(-3*t)*cos(9*exp(-3*t)*t)*(-6*exp(-3*t)*t+9*t^2*exp(-3*t))+1/2*exp(-3*t)*sin(9*exp(-3*t)*t)*(-6*exp(-3*t)*t+9*t^2*exp(-3*t))^2-1/2*exp(-3*t)*cos(9*exp(-3*t)*t)*(2*exp(-3*t)*t-12*t^2*exp(-3*t)+9*t^3*exp(-3*t))] 第二部分 数学问题求解与数据处理 1.(1) >> syms a t; >> f=sin(a*t)/t;laplace(f) ans = atan(a/s) (2) >> syms t a;f=t^5*sin(a*t);laplace(f) ans = 60*i*(-1/(s-i*a)^6+1/(s+i*a)^6) (3) >> syms t a; >> f=t^8*cos(a*t);laplace(f) ans = 20160/(s-i*a)^9+20160/(s+i*a)^9 2.(1) >> syms s a b;F=1/(s^2*(s^2-a^2)*(s+b));ilaplace(F) ans = 1/2/a^3/b^2/(a^2-b^2)*(2*t*a*b^3+2*(1-exp(-b*t)-b*t)*a^3+( -2*a+exp(a*t)*(a-b)+exp(-a*t)*(a+b))*b^2) (2) >> syms a b s;F=sqrt(s-a)-sqrt(s-b);ilaplace(F) ans = 1/2/t^(3/2)/pi^(1/2)*(exp(b*t)-exp(a*t)) (3) >> syms a b s; F=log((s-a)/(s-b));ilaplace(F) ans = (exp(b*t)-exp(a*t))/t 3 (1) >> syms x;f=x^2*(3*sym(pi)-2*abs(x));F=fourier(f) F = -6*(4+pi^2*dirac(2,w)*w^4)/w^4 >> ifourier(F) ans = x^2*(-4*x*heaviside(x)+3*pi+2*x) (2) >> syms t;f=t^2*(t-2*sym(pi))^2;F=fourier(f) F = 2*pi*(dirac(4,w)-4*pi^2*dirac(2,w)+4*i*pi*dirac(3,w)) >> ifourier(F) ans = x^2*(-2*pi+x)^2 4. (1) >> syms k a T;f=cos(k*a*T);F=ztrans(f) F = (z-cos(a*T))*z/(z^2-2*z*cos(a*T)+1) >> f1=iztrans(F) f1 = cos(a*T*n) (2) >> syms k T a;f=(k*T)^2*exp(-a*k*T);F=ztrans(f) F = T^2*z*exp(-a*T)*(z+exp(-a*T))/(z-exp(-a*T))^3 >> f1=iztrans(F) f1 = T^2*(1/exp(a*T))^n*n^2 (3) >> syms a k T;f=(a*k*T-1+exp(-a*k*T))/a; >> F=ztrans(f) F = 1/a*(a*T*z/(z-1)^2-z/(z-1)+z/exp(-a*T)/(z/exp(-a*T)-1)) >> iztrans(F) ans = (-1+a*T*n+(1/exp(a*T))^n)/a 5.用数值方法求解 (1) >> syms x; >> x1=solve('exp(-(x+1)^2+pi/2)*sin(5*x+2)') x1 = -2/5 验证过程 >> subs('exp(-(x+1)^2+pi/2)*sin(5*x+2)',x,x1) ans = 0 (2) >> syms x; >> y1=solve('(x^2+y^2+x*y)*exp(-x^2-y^2-x*y)=0','y') y1 = (-1/2+1/2*i*3^(1/2))*x (-1/2-1/2*i*3^(1/2))*x 验证过程 >> y2=simple(subs('(x^2+y^2+x*y)*exp(-x^2-y^2-x*y)=0','y',y1)) y2 = (x^2+(-1/2+1/2*i*3^(1/2))^2*x^2+x^2*(-1/2+1/2*i*3^(1/2)))*exp(-x^2-(-1/2+1/2*i*3^(1/2))^2*x^2-x^2*(-1/2+1/2*i*3^(1/2))) = 0 (x^2+(-1/2-1/2*i*3^(1/2))^2*x^2+x^2*(-1/2-1/2*i*3^(1/2)))*exp(-x^2-(-1/2-1/2*i*3^(1/2))^2*x^2-x^2*(-1/2-1/2*i*3^(1/2))) = 0 6. 首先求出积分: >> syms x c;y=int((exp(x)-c*x)^2,x,0,1) y = -1/2-2*c+1/2*exp(2)+1/3*c^2 编写一个出M文件: function y=new(c) y=-1/2-2*c+1/2*exp(2)+1/3*c^2; >> x=fminsearch('new',0) x = 3.0000 7. 编写M文件: function [c,ce]=f2(x); ce=[]; c=[x(1)+x(2);x(1)*x(2)-x(1)-x(2)+1.5;-10-x(1)*x(2)]; >> f=@(x)exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1); >> A=[];B=[];Aeq=[];Beq=[];xm=[-10;-10];xM=[10;10]; >> x0=(xm+xM)/2; >> ff=optimset;ff.TolX=1e-10;ff.TolFun=1e-20; >> x=fmincon(f,x0,A,B,Aeq,Beq,xm,xM,@f2,ff) Maximum number of function evaluations exceeded; increase OPTIONS.MaxFunEvals. x = 0.4195 0.4195 >> i=1;x=x0; >> while (1) [x,a,b]=fmincon('f',x,A,B,Aeq,Beq,xm,xM,'f2',ff); if b>0,break;end i=i+1; end >> x,i x = 1.1825 -1.7398 i = 5 8 书上所描述的ipslv_mex下载地址已经失效了，其他网站上也没有这个函数的下 载地址，所以这个题目没有找到该函数，运行失败。 >> f=-[592 381 273 55 48 37 23]; >> A=[3532 2356 1767 589 528 451 304];B=119567; >> intlist=[1;1;1;1;1;1;1];ctype=1; >> xm=zeros(7,1);xM=inf*ones(7,1); >> [res,b]=ipslv_mex(f,A,B,intlist,xM,xm,ctype) ??? Undefined command/function 'ipslv_mex'. . 9. >> syms x 通解的解法: >> y=dsolve('D2y-(2-1/x)*Dy+(1-1/x)*y=x^2*exp(-5*x)','x') y = exp(x)*C2+exp(x)*log(x)*C1+1/216*Ei(1,6*x)*exp(x)+11/1296*exp(-5*x)+5/216*exp(-5*x)*x+1/36*x^2*exp(-5*x) 特解的解法: >> syms x >> y=dsolve('D2y-(2-1/x)*Dy+(1-1/x)*y=x^2*exp(-5*x)', 'y(1)=syms(pi)','y(sym(pi))=1','x') y = -1/1296*exp(x)*(-1296*syms(pi)*exp(5)+6*exp(6)*Ei(1,6)+77)/exp(1)/exp(5)-1/1296*exp(x)*log(x)*(-1296*exp(1)*exp(5)+1296*exp(sym(pi))*syms(pi)*exp(5)-6*exp(sym(pi))*exp(6)*Ei(1,6)-77*exp(sym(pi))+6*exp(-5*sym(pi))*exp(6*sym(pi))*Ei(1,6*sym(pi))*exp(1)*exp(5)+11*exp(-5*sym(pi))*exp(1)*exp(5)+30*exp(-5*sym(pi))*sym(pi)*exp(1)*exp(5)+36*exp(-5*sym(pi))*sym(pi)^2*exp(1)*exp(5))/exp(sym(pi))/log(sym(pi))/exp(1)/exp(5)+1/1296*(6*exp(6*x)*Ei(1,6*x)+11+30*x+36*x^2)*exp(-5*x) 10. (1) >> syms t; >> x=dsolve('D2x+2*t*Dx+t^2*x=t+1') x = exp(t-1/2*t^2)*C2+exp(-t-1/2*t^2)*C1-1/2*i*pi^(1/2)*2^(1/2)* erf(1/2*i*2^(1/2)*(-1+t))*exp(-1/2+t-1/2*t^2) (2) >> syms x >> y=dsolve('Dy+2*x*y=x*exp(-x^2)','x') y = 1/2*(x^2+2*C1)*exp(-x^2) 11. >> f=@(t,x)[-x(2)-x(3);x(1)+0.2*x(2);0.2+(x(1)-5.7)*x(3)]; >> t_final=100; >> x0=[0;0;0]; >> [t,x]=ode45(f,[0,t_final],x0); 三维图 >>plot3(x(:,1),x(:,2),x(:,3));grid 25 20 15 10 5 0 10 15010 5-100 -58-20-10 6二维图: >> plot(x(:,1),x(:,2));grid 4 2 0 -2 -4 -6 -8 -10 -12 -10-5051015 改变: >> f=@(t,x)[-x(2)-x(3);x(1)+0.2*x(2);0.5+(x(1)-10)*x(3)]; >> t_final=100; >> x0=[0;0;0]; >> [t,x]=ode45(f,[0,t_final],x0); 三维图: >>plot3(x(:,1),x(:,2),x(:,3));grid 40 30 20 10 0 20 1020 100 0-10-10 -20-20 二维图: >> plot(x(:,1),x(:,2));grid 15 10 5 0 -5 -10 -15 -20 -15-10-505101520 12. 设x(1)=x,x(2)=dx,x(3)=y,x(4)=dy,x(5)=d2y >> f=inline(['[x(2);-x(1)-x(3)-(3*x(2))^2+(x(4))^3+6*x(5)+2*t;','x(4 );x(5);-x(5)-x(2)-exp(-x(1))-t]'],'t','x'); >> [t1,x1]=ode45(f,[1,0],[2,-4,-2,7,6]'); >> [t2,x2]=ode45(f,[1,2],[2,-4,-2,7,6]'); >> t=[t1(end:-1:1);t2];x=[x1(end:-1:1,:);x2]; >> plot(t,x) 12108 6 4 2 0 -2-4-68-800.20.40.60.811.21.41.61.82 6>>plot(x(:,1),x(:,3)) 4 2 0 -2-4-612345678910 13. >> [t,x,y]=sim('best',[0,10]);plot(t,x) 1.5 1 0.5 0 -0.5 -1012345678910 >> plot(t,y) 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8012345678910 14. >> t=0:0.2:2; >> y=t.^2.*exp(-5*t).*sin(t);plot(t,y,'o') 0.012 0.01 0.008 0.006 0.004 0.002 0t00.20.40.60.811.21.41.61.82x 10 >> ezplot('t.^2.*exp(-5*t).*sin(t)',[0,2]);hold on 10>> x1=0:0.01:2;y1=interp1(t,y,x1,'spline'); >> plot(x1,y1) 2-3 exp(-5 t) sin(t) 8 6 4 2 0 00.20.40.60.811.21.41.61.82 t