东北大学MATLAB实验参考答案.doc

东北大学MATLAB实验参考答案

《MATLAB语言与应用》实验课程任务书 信息科学与工程学院实验中心自动化实验室
《MATLAB语言与应用》实验课程任务书
一、 实验教学目标与基本要求
上机实验是本课程重要的实践教学环节；实验的目的不仅仅是验证理论知识，更重要的是通过上机实验，加强学生的实验手段与实践技能，掌握应用MATLAB语言求解问题的方法，培养学生分析问题、解决问题、应用知识的能力和创新精神，全面提高学生的综合素质。
上机实验共8学时。主要实验内容是基于理论课所学知识对课后典型习题进行MATLAB求解，基本掌握常见数学问题的求解方法与命令调用，更深入地认识和了解MATLAB语言强大的计算功能。
上机实验最终以书面报告的形式提交，并作为期末成绩考核内容的一部分。
二、 实验内容（8学时）
第一部分MATLAB语言编程、科学绘图与基本数学问题求解（4学时）
主要内容：掌握MATLAB语言编程基础、科学绘图方法、微积分问题、线性代数问题等基本数学问题的求解与应用。
练习题：
1、 安装MATLAB软件，应用demo命令了解主要功能，熟悉基本功能，会用help命令。
2、 用MATLAB语句输入矩阵和
，
前面给出的是矩阵，如果给出命令将得出什么结果？

Input
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];
A(5,6)=5
Answer=
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、 假设已知矩阵，试给出相应的MATLAB命令，将其全部偶数行提取出来，赋给矩阵，用命令生成矩阵，用上述命令检验一下结果是不是正确。

Input
A=magic(8);
B1=A(2:2:end, :)

Answer=
B1 =

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、 用数值方法可以求出，试不采用循环的形式求出和式的数值解。由于数值方法是采用double形式进行计算的，难以保证有效位数字，所以结果不一定精确。试采用运算的方法求该和式的精确值。

>> format long;sum(2.^[0:63])

ans =

1.844674407370955e+019

5、 选择合适的步距绘制出下面的图形。
（1），其中； （2），其中。
(1)>> t=-1:0.03:1; y=sin(1./t); plot(t,y)

>> t=[-1:0.03: -0.25, -0.248:0.001:0.248, 0.25:.03:1]; y=sin(1./t); plot(t,y)

(2)>> x=[-pi:0.05:pi];...
y=sin(tan(x))-tan(sin(x));...
plot(x,y)

x=[-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(x))-tan(sin(x));...
plot(x,y)

6、 试绘制出二元函数的三维图和三视图。

>> [x,y]=meshgrid(-2:.1:2);...
z=1./(sqrt((1-x).^2+y.^2))+1./(sqrt((1+x).^2+y.^2));...
surf(x,y,z),shading flat...

[x,y]=meshgrid(-2:.1:2);...
z=1./(sqrt((1-x).^2+y.^2))+1./(sqrt((1+x).^2+y.^2));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);

7、 试求出如下极限。
（1）； （2）； （3）。

（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);L=limit(limit(f,x,0),y,1)

L =

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、 已知参数方程，试求出和。

>> syms t; x=log(cos(t)); y=cos(t)-t*sin(t);
diff(y,t)/diff(x,t)

ans =

-(-2*sin(t)-t*cos(t))/sin(t)*cos(t)

>> f=diff(y,t,2)/diff(x,t,2); subs(f,t,sym(pi)/3)

ans =

3/8-1/24*pi*3^(1/2)

9、 假设，试求。
>> syms x y t
f=int(exp(-t^2),t,0,x*y);
x/y*diff(f,x,2)-2*diff(diff(f,x),y)+diff(f,y,2)
simple(ans)

ans =

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)

simplify:

-2*exp(-x^2*y^2)*(-x^2*y^2+1+x^3*y)

radsimp:

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)

combine(trig):

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)

factor:

-2*exp(-x^2*y^2)*(-x^2*y^2+1+x^3*y)

expand:

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)

combine:

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)

convert(exp):

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)

convert(sincos):

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)

convert(tan):

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)

collect(x):

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)

mwcos2sin:

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)

ans =

-2*exp(-x^2*y^2)*(-x^2*y^2+1+x^3*y)

10、 试求出下面的极限。
（1）；

>> syms k n; symsum(1/((2*k)^2-1),k,1,inf)

ans =

1/2
（2） 。
>> syms k n
limit(n*symsum(1/(n^2+k*pi),k,1,n),n,inf)

ans =

1
11、 试求出以下的曲线积分。
（1），为曲线，，
。

syms a t; x=a*(cos(t)+t*sin(t)); y=a*(sin(t)-t*cos(t));
f=x^2+y^2; I=int(f*sqrt(diff(x,t)^2+diff(y,t)^2),t,0,2*pi)

I =

2*csgn(a)*a^3*pi^2+4*csgn(a)*a^3*pi^4

（2），其中为正向上半椭圆。

>> syms x y a b c t; x=c*cos(t)/a; y=c*sin(t)/b;
P=y*x^3+exp(y); Q=x*y^3+x*exp(y)-2*y;
ds=[diff(x,t);diff(y,t)]; I=int([P Q]*ds,t,0,pi)

I =

-2/15*c*(-2*c^4+15*b^4)/b^4/a

12、 试求出Vandermonde矩阵的行列式，并以最简的形式显示 结果。

>> syms a b c d e; A=vander([a b c d e])
A =
[ 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(A), simple(ans)

ans =

(c-d)*(b-d)*(b-c)*(a-d)*(a-c)*(a-b)*(-d+e)*(e-c)*(e-b)*(e-a)

13、 试对矩阵进行Jordan变换，并得出变换矩阵。

>> 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、 试用数值方法和解析方法求取下面的Sylvester方程，并验证得出的结果。

15、 假设已知矩阵如下，试求出，，。

>> 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(-3*t)-1/2*t*exp(-3*t)+1/2*exp(-5*t)+1/2*t^2*exp(-3*t), 1/2*exp(-5*t)-1/2*exp(-3*t)+t*exp(-3*t), 1/2*t*exp(-3*t)+1/2*t^2*exp(-3*t), 1/2*exp(-5*t)-1/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(-5*t)-1/2*exp(-3*t), 1/2*exp(-3*t)+1/2*exp(-5*t), 1/2*t*exp(-3*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(-5*t)+1/2*exp(-3*t), -1/2*exp(-5*t)+1/2*exp(-3*t), exp(-3*t)+1/2*t*exp(-3*t), 1/2*t*exp(-3*t)-1/2*exp(-5*t)+1/2*exp(-3*t)]
[ -1/2*t^2*exp(-3*t), -t*exp(-3*t), -1/2*t^2*exp(-3*t)-t*exp(-3*t), exp(-3*t)-1/2*t^2*exp(-3*t)]

>> 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 x t; sin(A*t)

ans =

[ -sin(9/2*t), 0, sin(1/2*t), -sin(3/2*t)]
[ -sin(1/2*t), -sin(4*t), sin(1/2*t), -sin(1/2*t)]
[ sin(3/2*t), sin(t), -sin(5/2*t), sin(3/2*t)]
[ 0, -sin(t), -sin(t), -sin(3*t)]

>> 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 x t; exp(A*t)*sin(A^2*exp(A*t)*t)

ans =

[ exp(-9/2*t)*sin(t*(21*exp(-9/2*t)+2*exp(-1/2*t)-2*exp(3/2*t)+12))+sin(t*(5*exp(-9/2*t)+17*exp(-1/2*t)-3*exp(3/2*t)+5))+exp(1/2*t)*sin(t*(-11*exp(-9/2*t)-8*exp(-1/2*t)+6*exp(3/2*t)-11))+exp(-3/2*t)*sin(t*(-exp(-9/2*t)+6*exp(-1/2*t)+5*exp(3/2*t)+8)), exp(-9/2*t)*sin(t*(21+2*exp(-4*t)-2*exp(t)+12*exp(-t)))+sin(t*(5+17*exp(-4*t)-3*exp(t)+5*exp(-t)))+exp(1/2*t)*sin(t*(-11-8*exp(-4*t)+6*exp(t)-11*exp(-t)))+exp(-3/2*t)*sin(t*(-1+6*exp(-4*t)+5*exp(t)+8*exp(-t))), exp(-9/2*t)*sin(t*(23*exp(1/2*t)-2*exp(-5/2*t)+12*exp(-t)))+sin(t*(22*exp(1/2*t)-3*exp(-5/2*t)+5*exp(-t)))+exp(1/2*t)*sin(t*(-19*exp(1/2*t)+6*exp(-5/2*t)-11*exp(-t)))+exp(-3/2*t)*sin(t*(5*exp(1/2*t)+5*exp(-5/2*t)+8*exp(-t))), exp(-9/2*t)*sin(t*(21*exp(-3/2*t)+2*exp(-1/2*t)-2*exp(3/2*t)+12*exp(-3*t)))+sin(t*(5*exp(-3/2*t)+17*exp(-1/2*t)-3*exp(3/2*t)+5*exp(-3*t)))+exp(1/2*t)*sin(t*(-11*exp(-3/2*t)-8*exp(-1/2*t)+6*exp(3/2*t)-11*exp(-3*t)))+exp(-3/2*t)*sin(t*(-exp(-3/2*t)+6*exp(-1/2*t)+5*exp(3/2*t)+8*exp(-3*t)))]
[ exp(-1/2*t)*sin(t*(21*exp(-9/2*t)+2*exp(-1/2*t)-2*exp(3/2*t)+12))+exp(-4*t)*sin(t*(5*exp(-9/2*t)+17*exp(-1/2*t)-3*exp(3/2*t)+5))+exp(1/2*t)*sin(t*(-11*exp(-9/2*t)-8*exp(-1/2*t)+6*exp(3/2*t)-11))+exp(-1/2*t)*sin(t*(-exp(-9/2*t)+6*exp(-1/2*t)+5*exp(3/2*t)+8)), exp(-1/2*t)*sin(t*(21+2*exp(-4*t)-2*exp(t)+12*exp(-t)))+exp(-4*t)*sin(t*(5+17*exp(-4*t)-3*exp(t)+5*exp(-t)))+exp(1/2*t)*sin(t*(-11-8*exp(-4*t)+6*exp(t)-11*exp(-t)))+exp(-1/2*t)*sin(t*(-1+6*exp(-4*t)+5*exp(t)+8*exp(-t))), exp(-1/2*t)*sin(t*(23*exp(1/2*t)-2*exp(-5/2*t)+12*exp(-t)))+exp(-4*t)*sin(t*(22*exp(1/2*t)-3*exp(-5/2*t)+5*exp(-t)))+exp(1/2*t)*sin(t*(-19*exp(1/2*t)+6*exp(-5/2*t)-11*exp(-t)))+exp(-1/2*t)*sin(t*(5*exp(1/2*t)+5*exp(-5/2*t)+8*exp(-t))), exp(-1/2*t)*sin(t*(21*exp(-3/2*t)+2*exp(-1/2*t)-2*exp(3/2*t)+12*exp(-3*t)))+exp(-4*t)*sin(t*(5*exp(-3/2*t)+17*exp(-1/2*t)-3*exp(3/2*t)+5*exp(-3*t)))+exp(1/2*t)*sin(t*(-11*exp(-3/2*t)-8*exp(-1/2*t)+6*exp(3/2*t)-11*exp(-3*t)))+exp(-1/2*t)*sin(t*(-exp(-3/2*t)+6*exp(-1/2*t)+5*exp(3/2*t)+8*exp(-3*t)))]
[ exp(3/2*t)*sin(t*(21*exp(-9/2*t)+2*exp(-1/2*t)-2*exp(3/2*t)+12))+exp(t)*sin(t*(5*exp(-9/2*t)+17*exp(-1/2*t)-3*exp(3/2*t)+5))+exp(-5/2*t)*sin(t*(-11*exp(-9/2*t)-8*exp(-1/2*t)+6*exp(3/2*t)-11))+exp(3/2*t)*sin(t*(-exp(-9/2*t)+6*exp(-1/2*t)+5*exp(3/2*t)+8)), exp(3/2*t)*sin(t*(21+2*exp(-4*t)-2*exp(t)+12*exp(-t)))+exp(t)*sin(t*(5+17*exp(-4*t)-3*exp(t)+5*exp(-t)))+exp(-5/2*t)*sin(t*(-11-8*exp(-4*t)+6*exp(t)-11*exp(-t)))+exp(3/2*t)*sin(t*(-1+6*exp(-4*t)+5*exp(t)+8*exp(-t))), exp(3/2*t)*sin(t*(23*exp(1/2*t)-2*exp(-5/2*t)+12*exp(-t)))+exp(t)*sin(t*(22*exp(1/2*t)-3*exp(-5/2*t)+5*exp(-t)))+exp(-5/2*t)*sin(t*(-19*exp(1/2*t)+6*exp(-5/2*t)-11*exp(-t)))+exp(3/2*t)*sin(t*(5*exp(1/2*t)+5*exp(-5/2*t)+8*exp(-t))), exp(3/2*t)*sin(t*(21*exp(-3/2*t)+2*exp(-1/2*t)-2*exp(3/2*t)+12*exp(-3*t)))+exp(t)*sin(t*(5*exp(-3/2*t)+17*exp(-1/2*t)-3*exp(3/2*t)+5*exp(-3*t)))+exp(-5/2*t)*sin(t*(-11*exp(-3/2*t)-8*exp(-1/2*t)+6*exp(3/2*t)-11*exp(-3*t)))+exp(3/2*t)*sin(t*(-exp(-3/2*t)+6*exp(-1/2*t)+5*exp(3/2*t)+8*exp(-3*t)))]
[ sin(t*(21*exp(-9/2*t)+2*exp(-1/2*t)-2*exp(3/2*t)+12))+exp(-t)*sin(t*(5*exp(-9/2*t)+17*exp(-1/2*t)-3*exp(3/2*t)+5))+exp(-t)*sin(t*(-11*exp(-9/2*t)-8*exp(-1/2*t)+6*exp(3/2*t)-11))+exp(-3*t)*sin(t*(-exp(-9/2*t)+6*exp(-1/2*t)+5*exp(3/2*t)+8)), sin(t*(21+2*exp(-4*t)-2*exp(t)+12*exp(-t)))+exp(-t)*sin(t*(5+17*exp(-4*t)-3*exp(t)+5*exp(-t)))+exp(-t)*sin(t*(-11-8*exp(-4*t)+6*exp(t)-11*exp(-t)))+exp(-3*t)*sin(t*(-1+6*exp(-4*t)+5*exp(t)+8*exp(-t))), sin(t*(23*exp(1/2*t)-2*exp(-5/2*t)+12*exp(-t)))+exp(-t)*sin(t*(22*exp(1/2*t)-3*exp(-5/2*t)+5*exp(-t)))+exp(-t)*sin(t*(-19*exp(1/2*t)+6*exp(-5/2*t)-11*exp(-t)))+exp(-3*t)*sin(t*(5*exp(1/2*t)+5*exp(-5/2*t)+8*exp(-t))), sin(t*(21*exp(-3/2*t)+2*exp(-1/2*t)-2*exp(3/2*t)+12*exp(-3*t)))+exp(-t)*sin(t*(5*exp(-3/2*t)+17*exp(-1/2*t)-3*exp(3/2*t)+5*exp(-3*t)))+exp(-t)*sin(t*(-11*exp(-3/2*t)-8*exp(-1/2*t)+6*exp(3/2*t)-11*exp(-3*t)))+exp(-3*t)*sin(t*(-exp(-3/2*t)+6*exp(-1/2*t)+5*exp(3/2*t)+8*exp(-3*t)))]

第二部分 数学问题求解与数据处理（4学时）
主要内容：掌握代数方程与最优化问题、微分方程问题、数据处理问题的MATLAB求解方法。
练习题：
1、 对下列的函数进行Laplace变换。
（1） ；（2）；（3）。

(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、 对下面的式进行Laplace反变换。
（1） ；（2）；（3）。

(1)>> syms s a b; F=1/(s^2*(s^2-a^2)*(s+b)); ilaplace(F)

ans =

1/2/b^2/a^3/(a^2-b^2)*(2*t*a*b^3+2*(1-b*t-exp(-b*t))*a^3+(-2*a+exp(a*t)*(a-b)+(a+b)*exp(-a*t))*b^2)
(2)>> syms s a b; 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 =

1/t*(exp(b*t)-exp(a*t))

3、 试求出下面函数的Fourier变换，对得出的结果再进行Fourier反变换，观察是否能得出原来函数。
（1） ；（2）。

(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*(4*i*pi*dirac(3,w)-4*pi^2*dirac(2,w)+dirac(4,w))

>> ifourier(F)

ans =

x^2*(-2*pi+x)^2

4、 请将下述时域序列函数进行Z变换，并对结果进行反变换检验。
（1） ；（2）；（3）。

(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/exp(a*T))^n-1+a*T*n)/a
5、 用数值求解函数求解下述一元和二元方程的根，并对得出的结果进行检验。
（1） ；（2）。
(1) >> ezplot('exp(-(x+1)^2+pi/2)*sin(5*x+2)')

-2.93,-2.31,
(2)>> ezsurf('(x^2+y^2+x*y)*exp(-x^2-y^2-x*y)')

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

function y=exc6ff(c)
y=1/2*exp(1)^2+1/3*c^2-1/2-2*c;

>> x=fminsearch('exc6ff',0)

x =

3.00000000000000

7、 试求解下面的非线性规划问题。

function [c,ce]=exc6fun6a(x)
ce=[];
c=[x(1)+x(2); x(1)*x(2)-x(1)-x(2)+1.5; -10-x(1)*x(2)];

>> 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('exc6fun6',x0,A,B,Aeq,Beq,xm,xM,'exc6fun6a',ff)
Maximum number of function evaluations exceeded;
increase OPTIONS.MaxFunEvals

x =

0.41947326053910
0.41947326053910

8、 求解下面的整数线性规划问题。

function y=exc6fun2(x)
y=-(592*x(1)+381*x(2)+273*x(3)+55*x(4)+48*x(5)+37*x(6)+23*x(7));

>> f=[120 66 72 58 132 104];
A=[1 1 1 0 0 0; 0 0 0 1 1 1; 1 0 0 1 0 0; 0 1 0 0 1 0; 0 0 1 0 0 1];
B=[30; 18; 10; 18; 30]; intlist=[1;1;1;1;1];
ctype=[0;0;0;-1;1]; xm=zeros(5,1); xM=inf*ones(5,1);
[res,b]=ipslv_mex(f,A,B,intlist,xM,xm,ctype); res

res =

0
8
22
10
0
8

>> Aeq=[1 1 1 0 0 0; 0 0 0 1 1 1; 1 0 0 1 0 0]; Beq=[30; 18; 10];
A=[0 1 0 0 1 0; 0 0 -1 0 0 -1]; B=[18; -30];
intlist=ones(6,1); xm=zeros(6,1); xM=20000*ones(6,1); x0=xm;
[errmsg,f,x]=bnb20('exc6fun3',x0,intlist,xm,xM,A,B,Aeq,Beq);
if length(errmsg)==0, x=round(x), end

x =
0
8
22
10
0
8

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