4.2 微积分运算
4.2.1 导数
在MA TLAB 系统中为用户提供了一元显函数求导的符号计算函数diff ，可以调用此函数求符号导数，不但使用方便，而且计算准确，迅速，尤其是求结构复杂的高阶导数更显示出其优越性。用diff 可以求一元显函数的各阶导函数和在某点处的各阶导数。在这里将一元显函数)(x f y =的各阶导函数 ,''',",'y y y 分别记作 ,,,yxxx yxx yx 。用diff 作符号求导函数和在一点处的导数的调用格式和功能如表4-1。
表4-1 求一元显函数)(x f y =的各阶导函数的MA TLAB 方法
[例1] 设一元函数))25sin(23arctan(32-+=x x y ，求y 的一阶导函数)('x y 和四阶导函数)()
4(x y
。
解 下面用函数diff 的两种调用格式求解。
（1）输入计算y 对x 的一阶导函数和四阶导函数的程序：
>> syms x y
y=atan(3*x^2+2*sin(5*x^3-2)); yx=diff (y,x); yx4=diff (y,x,4);
y1=simple(yx), pretty(y1) y4=simple(yx4), pretty(y4)
运行后屏幕显示y 对x 的一阶导数和四阶导数如下：
y1 =
(6*x+30*cos(5*x^3-2)*x^2)/(1+(3*x^2+2*sin(5*x^3-2))^2) y4 =
(101250*sin(5*x^3-2)*x^8-81000*cos(5*x^3-2)*x^5-9000*s
in(5*x^3-2)*x^2)/(1+(3*x^2+2*sin(5*x^3-2))^2)-8*(-6750*cos(5*x ^3-2)*x^6-2700*sin(5*x^3-2)*x^3+60*cos(5*x^3-2))/(1+(3*x^2+2*s in(5*x^3-2))^2)^2*(3*x^2+2*sin(5*x^3-2))*(6*x+30*cos(5*x^3-2)*x^2)+48*(6-450*sin(5*x^3-2)*x^4+60*cos(5*x^3-2)*x)/(1+(3*x^2+2*sin(5*x^3-2))^2)^3*(3*x^2+2*sin(5*x^3-2))^2*(6*x+30*cos(5*x^3-2)*x^2)^2-12*(6-450*sin(5*x^3-2)*x^4+60*cos(5*x^3-2)*x)/(1+(3*x^2+2*sin(5*x^3-2))^2)^2*(6*x+30*cos(5*x^3-2)*x^2)^2-6*(6-450*sin(5*x^3-2)*x^4+60*cos(5*x^3-2)*x)^2/(1+(3*x^2+2*sin(5*x^3-2))^2)^2*(3*x^2+2*sin(5*x^3-2))-48*(6*x+30*cos(5*x^3-2)*x^2)^4/(1+(3*x^2+2*sin(5*x^3-2))^2)^4*(3*x^2+2*sin(5*x^3-2))^3+24*(6*x+30*cos(5*x^3-2)*x^2)^4/(1+(3*x^2+2*sin(5*x^3-2))^2)^3*(3*x^2+2*sin(5*x^3-2))
(2) 用下面的MATLAB 程序分别计算y 对x 的一阶导数和四阶导数：
>> syms x y
yx=diff((6*x+30*cos(5*x^3-2)*x^2)/(1+(3*x^2+2*sin(5*x^
3-2))^2),x)
yx4=diff((6*x+30*cos(5*x^3-2)*x^2)/(1+(3*x^2+2*sin(5*x
^3-2))^2),x,4)
所得到的结果与（1）相同。
[例2] 设一元函数6
51
2)613ln())25cos(23(2
5
32+-+
+-+=x x x x x y 。 （1）求y 的一阶导函数)('
x y 和五阶导函数)()5(x y ; （2）求y 在0=x 处的一阶导数)0('
y 和五阶导数)0()5(y 。 解 （1）输入计算y 对x 的一阶导函数)('
x y 和五阶导函数)()
5(x y
的程序：
>> syms x y
y=(3*x^2+2*cos(5*x^3-2))*log(13*x^5+6)+(2*x-1)/(5*x^2+
6);
yx=diff (y,x,1); yx5=diff (y,x,5); y1=simple(yx) y5=simple(yx5)
运行后屏幕显示化简后的y 对x 的一阶导函数y 1和五阶导函数y 5如下：
y1 =
(6*x-30*sin(5*x^3-2)*x^2)*log(13*x^5+6)+65*(3*x^2+2*co
s(5*x^3-2))*x^4/(13*x^5+6)+2/(5*x^2+6)-10*(2*x-1)/(5*x^2+6)^2*x
y5 =
2400000*(2*x-1)/(5*x^2+6)^5*x^3-507000*(6-450*cos(5*x^
3-2)*x^4-60*sin(5*x^3-2)*x)*x^7/(13*x^5+6)^2-2028000*(6*x-30*s in(5*x^3-2)*x^2)*x^6/(13*x^5+6)^2+325*(101250*cos(5*x^3-2)*x^8+81000*sin(5*x^3-2)*x^5-9000*cos(5*x^3-2)*x^2)*x^4/(13*x^5+6)+7800*(6-450*cos(5*x^3-2)*x^4-60*sin(5*x^3-2)*x)*x^2/(13*x^5+6)+7800*(6*x-30*sin(5*x^3-2)*x^2)*x/(13*x^5+6)-90000*(2*x-1)/(5*x^2+6)^4*x-360000/(5*x^2+6)^4*x^2+(-1518750*sin(5*x^3-2)*x^10+2025000*cos(5*x^3-2)*x^7+540000*sin(5*x^3-2)*x^4-18000*cos(5*x ^3-2)*x)*log(13*x^5+6)+2400000/(5*x^2+6)^5*x^4+6000/(5*x^2+6)^3+197730000*(3*x^2+2*cos(5*x^3-2))*x^10/(13*x^5+6)^3+2600*(6750*sin(5*x^3-2)*x^6-2700*cos(5*x^3-2)*x^3-60*sin(5*x^3-2))*x^3/(13*x^5+6)+1560*(3*x^2+2*cos(5*x^3-2))/(13*x^5+6)-42250*(6750*sin(5*x^3-2)*x^6-2700*cos(5*x^3-2)*x^3-60*sin(5*x^3-2))*x^8/(13*x^5+6)^2+5492500*(6-450*cos(5*x^3-2)*x^4-60*sin(5*x^3-2)*x)*x^12/(13*x^5+6)^3+65910000*(6*x-30*sin(5*x^3-2)*x^2)*x^11/(13*x^5+6)^3-535518750*(6*x-30*sin(5*x^3-2)*x^2)*x^16/(13*x^5+6)^4-4284150000*(3*x^2+2*cos(5*x^3-2))*x^15/(13*x^5+6)^4+27846975000*(3*x^2+2*cos(5*x^3-2))*x^20/(13*x^5+6)^5-12000000*(2*x-1)/(5*x^2+6)^6*x^5-2535000*(3*x^2+2*cos(5*x^3-2))*x^5/(13*x^5+6)^2
（2）输入计算y 在0 x 处的一阶导数)0('
y 和五阶导数)0()5(y 的程序：
>> x=0;
y1=(6*x-30*sin(5*x^3-2)*x^2)*log(13*x^5+6)+65*(3*x^2+2
*cos(5*x^3......-2))*x^4/(13*x^5+6)+2/(5*x^2+6)-10*(2*x-1)/(5*x^2+6)^2*x,
y5=2400000*(2*x-1)/(5*x^2+6)^5*x^3-507000*(6-450*cos(5
*x^3-2)*x^4......-60*sin(5*x^3-2)*x)*x^7/(13*x^5+6)^2-2028000*(6*x-30*sin(5*x^3-2)*x^2)*x^6/(13*x^5+6)^2+325*(101250*cos(5*x ^3-2)*x^8+81000*sin(5*x^3-2)*x^5-9000*cos(5*x^3-2)*x^2)*x^4/(13*x^5+6)+7800*(6-450*cos(5*x^3-2)*x^4-60*sin(5*x^3-2)*x)*x^2/(13*x^5+6)+7800*(6*x-30*sin(5*x^3-2)*x^2)*x/(13*x^5+6)-90000*(2*x-1)/(5*x^2+6)^4*x-360000/(5*x^2+6)^4*x^2+(-1518750*sin(5*x^3-2)*x^10+2025000*cos(5*x^3-2)*x^7+540000*sin(5*x^3-2)*x^4-18000*cos(5*x^3-2)*x)*log(13*x^5+6)+2400000/(5*x^2+6)^5*x^4+6000/(5*x^2+6)^3+197730000*(3*x^2+2*cos(5*x^3-2))*x^10/(13*x^5+6)^3+2600*(6750*sin(5*x^3-2)*x^6-2700*cos(5*x^3-2)*x^3-60*sin(5*x^3-2))*x^3/(13*x^5+6)+1560*(3*x^2+2*cos(5*x^3-2))/(13*x^5+6)-42250*(6750*sin(5*x^3-2)*x^6-2700*cos(5*x^3-2)*x^3-60*sin(5*x^3-2))*x^8/(13*x^5+6)^2+5492500*(6-450*cos(5*x^3-2)*x^4-60*sin(5*x ^3-2)*x)*x^12/(13*x^5+6)^3+65910000*(6*x-30*sin(5*x^3-2)*x^2)*x^11/(13*x^5+6)^3-535518750*(6*x-30*sin(5*x^3-2)*x^2)*x^16/(13*x^5+6)^4-4284150000*(3*x^2+2*cos(5*x^3-2))*x^15/(13*x^5+6)^4+27846975000*(3*x^2+2*cos(5*x^3-2))*x^20/(13*x^5+6)^5-12000000*(2*x-1)/(5*x^2+6)^6*x^5-2535000*(3*x^2+2*cos(5*x^3-2))*x^5/(13
*x^5+6)^2
运行后屏幕显示y 在0=x 处的一阶导数y 1和五阶导数y 5如下：
y1 =
0.3333 y5 =
-188.6186
[例3] 求下列函数的导数
(1) 2
53)
76(623213x x x x y +-?-= ; (2) 2
)5(ln +=x x y 。 解 输入程序：
>> syms x
y=(log(5*x))^(x+2);
f=3*x^3*(sqrt((23-6*x^5)/(6+7*x)^2))/(1-2*x); yx=diff (y,x); fx=diff(f,x); y1=simple(yx) f1=simple(fx)
运行后屏幕显示（1）和（2）的导数分别为f1和y1如下：
y1 =
log(5*x)^(x+2)*(log(log(5*x))+(x+2)/x/log(5*x)) f1 =
3*(294*x^7+135*x^6-198*x^5-322*x^2-230*x+414)*x^2/(
23-6*x^5)^(1/2)/(-1+2*x)^2/(6+7*x)^2
4.2.2 积分
如果函数)(x f 在区间],[b a 上连续，且)(x F 是)(x f 的一个原函数，则函数)(x f 在区间],[b a 上的定积分为
?
b
a
x f )(d )()(a F b F x -=。我们可以用MATLAB 软件中
的函数int 对定积分进行符号计算的调用格式和功能如表4-2所示。
[例4] 求
?
---5
4
)
3)(2)(1(5
t t t d t
解 输入程序：
>> syms t
f=5/((t-1)*(t-2)*(t-3)); F=int(f,t,4,5) y= double (F)
运行后屏幕显示计算的定积分值F 和近似值y 如下：
F =
25/2*log(2)-15/2*log(3) y =
0.4247
即 ?
---5
4
)
3)(2)(1(5
t t t d ≈t 0.4247。
[例5] 求
?
++1
0n
e
dx b ax d x
解 输入程序：
>> syms x a b d e n
f=(a*x+b)/((d*x+e)^(1/n)); F=int(f,0,1) y=simple(F)
运行后屏幕显示计算的定积分F 和化简后的定积分y 如下：
F =
n*(-2*e*(d+e)^(1/n)*n*b*d+e^2*(d+e)^(1/n)*n*a+2*e^(1/n)
*n*b*e*d+e^(1/n)*a*d^2*n+2*e^(1/n)*n*b*d^2-e^(1/n)*n*a*e^2-e^(1/n)*b*d^2-e^(1/n)*d*a*e-e^(1/n)*a*d^2+e*(d+e)^(1/n)*b*d-e^(1/n)*e*b*d)/(e^(1/n))/((d+e)^(1/n))/d^2/(-3*n+1+2*n^2)
y =
n*((-2*e*(d+e)^(1/n)*b*d+e^2*(d+e)^(1/n)*a+e^(1/n)*a*d^
2+2*e^(1/n)*b*d^2-e^(1/n)*a*e^2+2*e^(1/n)*e*b*d)*n-e^(1/n)*a*d ^2+e*(d+e)^(1/n)*b*d-e^(1/n)*b*d^2-e^(1/n)*d*a*e-e^(1/n)*e*b*d )/(e^(1/n))/((d+e)^(1/n))/d^2/(-3*n+1+2*n^2)
[例6] 求下列不定积分
（1）
?
-2
1
2
292x
x d x ；（2）
?
-2
1
221x x d x ；（3）?+2132)
4(1
x d x 。
解 下面将三个不定积分的被积函数f1，f2，f3构成一个矩阵f ，用一个积分命
令求不定积分，即输入程序：
>> syms x
f1=(2*x^2)/(sqrt(9-x^2)); f2=sqrt(x^2-1) /(2*x); f3=1/(sqrt((4+x^2)^3)); f=[ f1; f2; f3]; F=int(f,1,2) y= double(F)
运行后屏幕显示三个定积分值的矩阵F 和近似值矩阵y 如下：
F =
[ -2*5^(1/2)+9*asin(2/3)+2*2^(1/2)-9*asin(1/3)] [ 1/2*3^(1/2)-1/6*pi] [ 1/8*2^(1/2)-1/20*5^(1/2)] y =
1.8653 0.3424 0.0650 即
?
-2
1
2
292x
x d ≈x 1.8653 ，
?
-2
1
221x x d ≈x 0.3424，?+2132)
4(1
x d ≈x 0.0650。
[例7] 求下列定积分 （1）
?
ππ+---+52322)1)(1(12t t t t t d t ； （2）?ππ
++523
2
12cos 1cos t t d t 。 解 下面将两个不定积分的被积函数f1，f2构成一个矩阵f ，用一个积分命令求
不定积分，即输入程序：
>> syms t
f1=(t^2+2*t-1)/((t-1)*(t^2-t+1)); f2=((cos(t))^2+1)/(cos(2*t)+1); f=[ f1; f2];
F=int(f,pi/3,2*pi/5) y= double(F)
运行后屏幕显示两个定积分的矩阵F 及其近似值的矩阵y 如下：
F =
[ -log(5)+2*log(2*pi-5)-1/2*log(25+4*pi^2-10*pi)+5/3*3^(
1/2)*atan((4/15*pi-1/3)*3^(1/2))+log(3)-2*log(pi-3)+1/2*log(9+pi^2-3*pi)-5/3*3^(1/2)*atan((2/9*pi-1/3)*3^(1/2))]
[1/30*(pi*5^(1/2)-pi+15*2^(1/2)*(5+5^(1/2))^(1/2)-15*3^(
1/2)*5^(1/2)+15*3^(1/2))/(5^(1/2)-1)]
y =
3.7172
0.7775
即
?
ππ+---+523
22)
1)(1(1
2t t t t t d t ≈3.7172，
?
ππ++523
21
2cos 1
cos t t d t ≈0.7775。
[例8] 计算下列定积分： （1）
?
1
e
t
2d t ；（2）?1
02
)3(cos t d t ;（3）?+1
02
22
)1(t t d t .
解 因为题设给定的三个定积分的上、下限分别相同，所以可以用向量的形式输入
程序：
>> syms t
A = [exp(sqrt(2*t)),(cos(3*t))^2,(t^2/((1+t^2)^2))]; g=int(A,t, 0,1) S=double(g)
运行后屏幕显示计算三个定积分的值的向量g 及其近似值的向量S 依次为
g =
[exp(2^(1/2))*2^(1/2)-exp(2^(1/2))+1,
1/6*cos(3)*sin(3)+1/2, 22894441/3063060-19/8*pi]
S =
2.7038 0.4767 0.0131 即
?
1
e t
2d ≈t 2.7038, ?1
02
)3(cos t d ≈t 0.4767和?+1
02
22
)1(t t d ≈t 0.0131。
[例9] 已知???
??>-+≤+=,1,752
1,1,sin 1)(2x x x x x x f 求?20)(x f d x 。
解 输入程序：
>> syms x
f1=sin(x)+1; f2=(x^2)/2+5*x-7; F=int(f1,x,0,1)+ int(f2,x,1,2) y= double (F)
运行后屏幕显示计算的定积分值F 和近似值y 如下：
F =
-cos(1)+11/3 y =
3.1264 即
?
2
)(x f d ≈x 3.1264。
[例10] 计算
()?π+0
cos 32x d x .
解 由于 ???
????
π≤<π-π≤≤=x x x x x 2,cos 32
0,cos 3cos 3
根据积分的可加性，有
()?π
+0
c o s
32x d ()?π
+=20
cos 32x x d ()?π
π-+2
cos 32x x d x . 输入程序：
>> syms x
f1=3*cos(x)+2; f2=-3*cos(x)+2;
F=int(f1,x,0,pi/2)+ int(f2,x,pi/2,pi) y= double (F)
运行后屏幕显示计算的定积分值F 和近似值y 如下：
F =
6+2*pi y =
12.2832
Matlab GUI编程教程（适用于初学者） 1．首先我们新建一个GUI文件：如下图所示； 选择Blank GUI(Default) 2．进入GUI开发环境以后添加两个编辑文本框，6个静态文本框，和一个按钮，布置如下
图所示； 布置好各控件以后，我们就可以来为这些控件编写程序来实现两数相加的功能了。3．我们先为数据1文本框添加代码； 点击上图所示红色方框，选择edit1_Callback，光标便立刻移到下面这段代码的位置。 1. 2. 3.function edit1_Callback(hObject, eventdata, handles) 4.% hObject handle to edit1 (see GCBO) 5.% eventdata reserved - to be defined in a future version of MATLAB
6.% handles structure with handles and user data (see GUIDATA) 7.% Hints: get(hObject,'String') returns contents of edit1 as text 8.% str2double(get(hObject,'String')) returns contents of edit1 as a double 复制代码 然后在上面这段代码的下面插入如下代码： 1. 2.%以字符串的形式来存储数据文本框1的内容. 如果字符串不是数字，则现实空白内容input = str2num(get(hObject,'String')); %检查输入是否为空. 如果为空,则默认显示为0if (isempty(input)) set(hObject,'String','0')endguidata(hObject, handles); 复制代码 这段代码使得输入被严格限制，我们不能试图输入一个非数字。 4．为edit2_Callback添加同样一段代码 5 现在我们为计算按钮添加代码来实现把数据1和数据2相加的目的。 用3中同样的方法在m文件中找到pushbutton1_Callback代码段 如下； 1.function pushbutton1_Callback(hObject, eventdata, handles) 2.% hObject handle to pushbutton1 (see GCBO) 3.% eventdata reserved - to be defined in a future version of MATLAB 4.% handles structure with handles and user data (see GUIDATA) 复制代码
实例matlab-非线性规划-作业
现代设计方法-工程优化理论、方法与设计 姓名 学号 班级 研 问题 ： 某厂向用户提供发动机，合同规定，第一、二、三季度末分别交货40台、60台、80台。每季度的生产费用为 （元），其中x 是该季生产的台数。若交货后有剩余，可用于下季度交货，但需支付存储费，每台每季度c 元。已知工厂每季度最大生产能力为100台，第一季度开始时无存货，设a=50、b=0.2、c=4，问工厂应如何安排生产计划，才能既满足合同又使总费用最低。讨论a 、b 、c 变化对计划的影响，并作出合理的解释。 问题的分析和假设： 问题分析：本题是一个有约束条件的二次规划问题。决策变量是工厂每季度生产的台数，目标函数是总费用（包括生产费用和存储费）。约束条件是生产合同，生产能力的限制。在这些条件下需要如何安排生产计划，才能既满足合同又使总费用最低。 问题假设： 1、工厂最大生产能力不会发生变化； 2、合同不会发生变更； 3、第一季度开始时工厂无存货； 4、生产总量达到180台时，不在进行生产； 5、工厂生产处的发动机质量有保证，不考虑退货等因素； 6、不考虑产品运输费用是否有厂家承担等和生产无关的因素。 符号规定： x1——第一季度生产的台数； x2——第二季度生产的台数； 180-x1-x2——第三季度生产的台数； y1——第一季度总费用； y2——第二季度总费用； y3——第三季度总费用； y ——总费用（包括生产费用和存储费）。 ()2bx ax x f +=
建模： 1、第一、二、三季度末分别交货40台、60台、80台； 2、每季度的生产费用为 （元）； 3、每季度生产数量满足40 ≤x1≤100，0≤x2≤100，100≤x1+x2 ≤180； 4、要求总费用最低，这是一个目标规划模型。 目标函数: y1 2111x b x a Z ?+?= y2()4012222-?+?+?=x c x b x a Z y3()()()10018018021221213 -+?+--?+--?=x x c x x b x x a Z y x x x x x x Z Z Z Z 68644.04.04.0149201 212221321--+++=++= 40≤x1≤100 0≤x2≤100 100≤x1+x2≤180 ()2 bx ax x f +=
MATLAB简单程序大全 求特征值特征向量 A=[2 3 4;1 5 9;8 5 2] det(A) A' rank(A) inv(A) rref(A) eig(A)%求特征值和特征向量 卫星运行问题 h=200,H=51000,R=6378; a=(h+H+2*R)/2; c=(H-h)/2; b=(a^2-c^2)^(1/2); e=c/a; f=sqrt(1-exp(2).*cos(t)^2); l=int(f,t,0,pi/2) L=4*a.*l 动态玫瑰线 n=3;N=10000; theta=2*pi*(0:N)/N; r=cos(n*theta); x=r.*cos(theta); y=r.*sin(theta); comet(x,y) 二重积分 syms x y f=x^2*sin(y); int(int(f,x,0,1),y,0,pi) ezmesh(f,[0,1,0,pi]) 函数画图 syms x;f=exp(-0.2*x)*sin(0.5*x); ezplot(f,[0,8*pi])
玫瑰线 theta=0:0.01:2*pi; r=cos(3*theta); polar(theta,r,'r') 求x^2+y^2=1和x^2+z^2=1所围成的体积 syms x y z R r=1; Z=sqrt(1-x^2); y0=Z; V=8*int(int(Z,y,0,y0),x,0,1) 求导数及图像 f='1/(5+4*cos(x))'; subplot(1,2,1);ezplot(f) f1=diff(f) subplot(1,2,2);ezplot(f1) 绕x轴旋转 t=(0:20)*pi/10; r=exp(-.2*t).*sin(.5*t); theta=t; x=t'*ones(size(t)); y=r'*cos(theta); z=r'*sin(theta); mesh(x,y,z) colormap([0 0 0]) 某年是否闰年 year=input('input year:='); n1=year/4; n2=year/100; n3=year/400; if n1==fix(n1)&n2~=fix(n2) disp('是闰年') elseif n1==fix(n1)&n3==fix(n3) disp('是闰年') else
p=p1';t=t1'; [pn,minp,maxp,tn,mint,maxt]=premnmx(p,t); %原始数据归一化 net=newff(minmax(pn),[5,1],{'tansig','purelin'},'traingdx'); %设置网络,建立相应的BP网络 net.trainParam.show=2000; % 训练网络 net.trainParam.lr=0.01; net.trainParam.epochs=100000; net.trainParam.goal=1e-5; [net,tr]=train(net ,pn,tn); %调用TRAINGDM算法训练BP网络 pnew=pnew1'; pnewn=tramnmx(pnew,minp,maxp); anewn=sim(net,pnewn); %对BP网络进行仿真 anew=postmnmx(anewn,mint,maxt); %还原数据 y=anew'; 1、BP网络构建 （1）生成BP网络 = net newff PR S S SNl TF TF TFNl BTF BLF PF (,[1 2...],{ 1 2...},,,) R?维矩阵。 PR：由R维的输入样本最小最大值构成的2 S S SNl：各层的神经元个数。 [ 1 2...] { 1 2...} TF TF TFNl：各层的神经元传递函数。 BTF：训练用函数的名称。 （2）网络训练 [,,,,,] (,,,,,,) = net tr Y E Pf Af train net P T Pi Ai VV TV （3）网络仿真 = [,,,,] (,,,,) Y Pf Af E perf sim net P Pi Ai T {'tansig','purelin'},'trainrp'
Matlab 程序代码 绘制y = 10e-1.5t sin( 7.75t ) 的函数图象 7.75 fv clear; t=0:0.02:10; f1=10/sqrt(7.75).*exp(-1.5*t); f2=sin(sqrt(7.75).*t); y=f1.*f2; plot(t,y,'-k',t,y,'ok'); xlabel('t');ylabel('y(t) ');title('函数图像') axis([-2 10 -0.5 2]) 拉氏变换 clear; clc; syms s t fs1 fs2 fs3 ft1 ft2 ft3; L=1,C=0.1,R=[1.5 3 5]; h1=1/(L*C*s^2+R(1)*C*s+1); h2=1/(L*C*s^2+R(2)*C*s+1); h3=1/(L*C*s^2+R(3)*C*s+1); fs1=h1*(1/s); fs2=h2*(1/s); fs3=h3*(1/s);
ft1=ilaplace(fs1,s,t); ft2=ilaplace(fs2,s,t); ft3=ilaplace(fs3,s,t); ezplot(t,ft1); hold on; ezplot(t,ft2); hold on; ezplot(t,ft3); 信号编码 对[1 1 0 1 1 1 0 1 0 0 1]进行编码。 clear; clc; c=[1 1 0 1 1 1 0 1 0 0 1] for i=1:length(c) if i==1 d1(i)=0;d2(i)=0; elseif i==2 d1(i)=c(i-1);d2(i)=c(i-1); elseif i==3 d1(i)=mod(c(i-1)+c(i-2),2); d2(i)=c(i-1); else d1(i)=mod(c(i-1)+c(i-2),2); d2(i)=mod(c(i-1)+c(i-3),2); end
MATLAB作业一1、试求出如下极限。 （1） 23 25 (2)(3) lim (5) x x x x x x x ++ + →∞ ++ + ，（2） 23 3 1 2 lim () x y x y xy x y →- → + + ，（3） 22 22 22 1cos() lim ()x y x y x y x y e+ → → -+ + 解：（1）syms x； f=((x+2)^(x+2))*((x+3)^(x+3))/((x+5)^(2*x+5)) limit（f，x，inf） =exp（-5） （2）syms x y； f=(x^2*y+x*y^3)/(x+y)^3； limit(limit(f,x,-1),y,2) =-6； （3）syms x y； f=(1-cos(x^2+y^2))/(x^2+y^2)*exp(x^2+y^2)； limit(limit(f,x,0),y,0) =0 2、试求出下面函数的导数。 （1 ）() y x=, (2)22 atan ln() y x y x =+ 解; （1）syms x; f=sqrt(x*sin(x)*sqrt(1-exp(x))); g= diff(f,x); g== (sin(x)*(1 - exp(x))^(1/2) + x*cos(x)*(1 - exp(x))^(1/2) - (x*exp(x)*sin(x))/(2*(1 - exp(x))^(1/2)))/(2*(x*sin(x)*(1 - exp(x))^(1/2))^(1/2)) pretty(g)= (2)syms x y; f=atan(y/x)-log(x^2+y^2) pretty(-simple(diff(f,x)/diff(f,y)))= 2 x + y =------- x - 2 y (3) 假设1 cos u- =，试验证 22 u u x y y x ?? = ???? 。 解：syms x y; u=1/cos(sqrt(x/y)); diff(diff(u,x),y)-diff(diff(u,y),x)=0; 所以： 22 u u x y y x ?? = ????
matlab 仿真实例 实验五MATLAB 及仿真实验一、控制系统的时域分析 (一)稳定性 1、系统传递函数为G(s)，试判断其稳定性。 程序： >> nu m=[3,2,5,4,6]; >> den=[1,3,4,2,7,2]; >> sys=tf( nu m,de n); >> figure(1); >> pzmap(sys); >> title(' 零极点图') 由图可知：在S 右半平面有极点，因此可知系统是不稳定的。 2、用MATLA 求 出 G(s)=(s A 2+2*s+2)/(s A 4+7*s A 3+5*s+2) 的极点。 程序及结果： >> sys=tf([1,2,2],[1,7,3,5,2]); >> p=pole(sys) 矿'. 赳 _ ■ —
-6.6553 0.0327 + 0.8555i 0.0327 - 0.8555i -0.4100 (二)阶跃响应 1、二阶系统G(s)=10/s A2+2*s+10 1)键入程序，观察并记录单位阶跃响应曲线： 程序： >> sys=tf(10,[1,2,10]); >> step(sys); >> title('G(s)=10/sA2+2*s+10 单位阶跃响应曲线') 2）计算系统闭环跟、阻尼比、无阻尼振荡频率，并记录程序及结果： >> sys=tf(10,[1,2,10]); >> p=pole(sys)
p = -1.0000 + 3.0000i -1.0000 - 3.0000i >> [wn,z]=damp(sys) wn = 3.1623 3.1623 z = 0.3162 0.3162 3)记录实际测取的峰值大小，峰值时间和过渡过程时间，并填表实际值理论值峰值Cmax 1.35s 峰值时间tp 1.05s 过渡时间+5% 3.54s ts +2% 3.18s 程序: >> sys=tf(10,[1,2,10]); >> step(sys); >> title('G(s)=10/sA2+2*s+10 单位阶跃响应曲线')
MATLAB 作业二 1、请将下面给出的矩阵A 和B 输入到MATLAB 环境中，并将它们转换成符号矩阵。若某一 矩阵为数值矩阵，另以矩阵为符号矩阵，两矩阵相乘是符号矩阵还是数值矩阵。 57651653 5501232310014325462564206441211346,3 9636623 51521210760077410120172440773 473 78 124867217110 7 681 5A B ???? ?????????? ????? ?==?????? ????? ?---????????--??? ? 解：A 转换为符号矩阵；a=sym（A） a=[5,7,6,5,1,6,5] [2,3,1,0,0,1,4][6,4,2,0,6,4,4][3,9,6,3,6,6,2][10,7,6,0,0,7,7][7,2,4,4,0,7,0][4,8,6,7,2,1,7]B 转换为符号矩阵；b=sym（B）b = [3,5,5,0,1,2,3][3,2,5,4,6,2,5][1,2,1,1,3,4,6][3,5,1,5,2,1,2][4,1,0,1,2,0,1][-3,-4,-7,3,7,8,12][1,-10,7,-6,8,1,5] 若某一矩阵为数值矩阵，另以矩阵为符号矩阵，两矩阵相乘是符号矩阵例；a*B= [48,3,64,48,159,106,194][17,-26,47,-8,62,26,59][48,-8,52,12,108,64,124][59,22,41,69,151,101,184][43,-22,91,13,175,121,220][22,39,4,53,88,94,147][75,11,115,36,151,70,151] 2、利用MATLAB 语言提供的现成函数对习题1中给出的两个矩阵进行分析，判定它们是否 为奇异矩阵，得出矩阵的秩、行列式、迹和逆矩阵，检验得出的逆矩阵是否正确。 解：由于a=det（A）=3.7396e+04；故A 是非奇异矩阵。B=det（B）=0，故B 是奇异矩阵； 由于a=rank(A)=7，故A 的秩为7；由于b=rank(B)=5，故B 的秩为5；由于a=trace(A)=27，b=trace(B)=26，故A，B 的迹为27,26；由a=inv(A)得A 的逆矩阵如下；
Matlab 程序代码 绘制 1.5 10sin(7.75)7.75 t y e t -= 的函数图象。 fv clear; t=0:0.02:10; f1=10/sqrt(7.75).*exp(-1.5*t); f2=sin(sqrt(7.75).*t); y=f1.*f2; plot(t,y,'-k',t,y,'ok'); xlabel('t');ylabel('y(t) ');title('函数图像') axis([-2 10 -0.5 2]) 拉氏变换 clear; clc; syms s t fs1 fs2 fs3 ft1 ft2 ft3; L=1,C=0.1,R=[1.5 3 5]; h1=1/(L*C*s^2+R(1)*C*s+1); h2=1/(L*C*s^2+R(2)*C*s+1);
h3=1/(L*C*s^2+R(3)*C*s+1); fs1=h1*(1/s); fs2=h2*(1/s); fs3=h3*(1/s); ft1=ilaplace(fs1,s,t); ft2=ilaplace(fs2,s,t); ft3=ilaplace(fs3,s,t); ezplot(t,ft1); hold on; ezplot(t,ft2); hold on; ezplot(t,ft3); 信号编码 对[1 1 0 1 1 1 0 1 0 0 1]进行编码。 clear; clc; c=[1 1 0 1 1 1 0 1 0 0 1] for i=1:length(c) if i==1 d1(i)=0;d2(i)=0; elseif i==2 d1(i)=c(i-1);d2(i)=c(i-1); elseif i==3
《Matlab/Simulink电力系统建模与仿真》 上机实验报告 班级：15电气工程及其自动化二班 学号：154139240096 姓名：汤嘉旺 实验一：Powergui在简单电力系统潮流计算中的应用实例 1.1实验内容与要求 完成2机5节点电力系统的潮流计算，以2机5节点电力系统为模型进行Matlab/Simulink电力系统建模与仿真。并完成电力系统元件的模型选择、模型参数的计算及设置、计算结果及比较。 2机5节点电力系统图 1.2 Simulink建模原理图和主要模块参数设置 1、发电机模型 在该系统中的两台发电机均选用p.u.标准同步电机模块
“Synchronous Machine pu Standard”,该模块使用标幺值参数，以转子dq 轴建立的坐标系为参数，定子绕组为星形连接。 2、变压器模块 系统中的两台变压器均选用三相双绕组变压器模块 “Three-phase Transformer(Two Windings)”,采用Y-Y连接方式。 3、线路模块 系统中带有地导纳的线路选用三相“II”形等值模块 “Three Phase PI Section Line”,没有对地导纳的线路选用三相串联RLC支路模块“Three Phase Series RLC Branch”。 4、负荷模块 在SimPowerSystems库中，利用R、L、C的串联或并联组合，提供了两个静态三相负荷模块。这两种模块是用恒阻抗支路模拟负荷，在仿真时，在给定的频率下负荷阻抗为常数。 5、母线模块 选择带有测量元件的母线模型，及三相电压电流测量元 件“Three-Phase V-I Measurement”来模拟系统中的母线，同时方便测量流过线路的潮流，在线路元件的两端也设置了该元件。 2机5节点电力系统潮流计算仿真模型图
Matlab 程序代码 绘 制 1.5sin(7.75)7.75t y e t -=的函数图象。 fv clear; t=0:0.02:10; f1=10/sqrt(7.75).*exp(-1.5*t); f2=sin(sqrt(7.75).*t); y=f1.*f2; plot(t,y,'-k',t,y,'ok'); xlabel('t');ylabel('y(t) ');title('函数图像') axis([-2 10 -0.5 2]) 拉氏变换 clear; clc; syms s t fs1 fs2 fs3 ft1 ft2 ft3; L=1,C=0.1,R=[1.5 3 5]; h1=1/(L*C*s^2+R(1)*C*s+1); h2=1/(L*C*s^2+R(2)*C*s+1);
h3=1/(L*C*s^2+R(3)*C*s+1); fs1=h1*(1/s); fs2=h2*(1/s); fs3=h3*(1/s); ft1=ilaplace(fs1,s,t); ft2=ilaplace(fs2,s,t); ft3=ilaplace(fs3,s,t); ezplot(t,ft1); hold on; ezplot(t,ft2); hold on; ezplot(t,ft3); 信号编码 对[1 1 0 1 1 1 0 1 0 0 1]进行编码。clear; clc; c=[1 1 0 1 1 1 0 1 0 0 1] for i=1:length(c) if i==1 d1(i)=0;d2(i)=0; elseif i==2 d1(i)=c(i-1);d2(i)=c(i-1); elseif i==3
Matlab 综合作业：运用matlab 建模解题 导弹追踪问题 1. 位于坐标原点的甲舰向位于x 轴上点A (1, 0)处的乙舰发射导弹， 导弹头始终对准乙舰．如果乙舰以最大的速度v 0(常数)沿平行于y 轴的直线行驶，导弹的速度是5v 0，求导弹运行的曲线方程．乙舰行驶多远时，导弹将它击中？ 解：假设t 时刻导弹的位置为P (x (t ), y (t ))，乙舰位于),1(0t v Q 由于导弹头 始终对准乙舰，故此时直线PQ 就是导弹的轨迹曲线弧OP 在点P 处的切线， 即有 x y t v y --=1'0 即 y y x t v +-=')1(0 (1) 又根据题意，弧OP 的长度为AQ 的5倍， 即 0d 5x x v t =? (2) 由(1),(2)消去t, 整理得模型: (3) '151 ")1(2y y x +=- 值条件为: 0)0(=y 0)0('=y 令y 1=y , y 2=y 1`，将方程（3）化为一阶微分方程组． 2 151 '')1(y y x +=- ? ?????-+==)1/(151''21221x y y y y 1.建立M 文件a1.m function dy=a1(x,y) dy=zeros(2,1); dy(1)=y(2); dy(2)=1/5*sqrt(1+y(1)^2)/(1-x); 2. 取x0=0，xf=0．9999，建立主程序ff6．m 如下： x0=0,xf=0.9999 [x,y]=ode15s('a1',[x0 xf],[0 0]);
plot(x,y(:,1),'b.') hold on y=0:0.01:2; plot(1,y,'b*') 运行得图： 结论: 导弹大致在（1，0．2）处击中乙舰.
最优化平时作业 一、目标规划 1、题目：见书中例题P110 例 4 2、解题方法：利用Lingo 求解 3、具体步骤 （1）.对应于第一优先等级，建立线性规划问题： model: min=-d1; 5*x1+10*x2<=60; x1-2*x2+d1_-d1=0; end 运行结果：-d1=0 （2 ）对应于第二优先等级，将-d1=0 作为约束条件，建立线性规划问题：min=d2_; 5*x1+10*x2<=60; x1-2*x2+d1_-d1=0; 4*x1+4*x2+d2_-d2=36; -d1=0; end 运行结果：d2=0; （3）.对应于第三优先等级，将-d1=0, -d1=0 作为约束条件，建立线性规划问题：min=d3_; 5*x1+10*x2<=60; x1-2*x2+d1_-d1=0; 4*x1+4*x2+d2_-d2=36; 6x1+8*x2+d3_-d3=48; -d1=0; d2=0; end 运行结果：d3=0; X1 4.800000 X2 2.400000 二、动态规划之0-1 背包问题 1、题目：给定n种物品和一背包。物品i 的重量是Wi ，其价值为Vi ，背包的容量是c，问应如何选择装入背包中的物品，使得装入背包中物品的总价值最大。 2、解题方法与思路：利用java 求解，.思想方法是回溯思想 3、需求分析 对于给定n 种物品和一背包。在容量最大值固定的情况下，要求装入的物品价值最大化 4、java 源程序及运行结果 BackTrace.java * To change this template, choose Tools | Template Manager * and open the template in the editor. */ package sunfa; import java.util.Date;
一个实例搞定MATLAB界面编程 作者：彭军 邮件：pjun9@https://www.doczj.com/doc/e34773373.html, 博客：https://www.doczj.com/doc/e34773373.html,/pengjun 下面请跟我一步一步做一个图像处理的程序，如果您坚持做完这个实例，我想MATLAB界面编程对您而言，就没有什么难度了。当然，我这里说的是，您首先要有一定的MATLAB 编程基础。还有，我的MATLAB版本是2008a。在2008a以前的版本中没有工具栏编辑器，如果需要工具栏要手动写程序，这个我就不多讲了。好了，废话少说，跟我来吧！ 在MATLAB的命令窗口(Command Window)中运行guide命令，来打开GUIDE界面，如下： 然后，选择空模板(Blang GUI)，点击OK，即可打开GUIDE的设计界面，如下：
点击工具栏上的菜单编辑器(Menu Editor)，打开菜单编辑器，如下： 在Menu Bar中新建一个菜单项，名字为“文件”，其他设置请看下图： 在“文件”菜单下添加菜单项：“打开”，“保存”，“退出”。见下图：
如果需要在菜单项“退出”上面添加一个分割线的话，选中“Separator above this item”就行了。
保存我的界面为pjimage.fig.保存完毕之后，会自动打开pjimage.m 文件，而我们所有的程序都是要写在这个M 文件里面的。在编程中，我们的每一个鼠标动作都对应一个Callback 函数。那么我们的菜单项也是如此的。 在界面上，单击鼠标右键选择“Property Inspector ”，即可打开属性窗口。当我们点击不同的 然后，点击工具栏的保存按钮。之后，点击工具栏的运行按钮(Run Figure)。注意，工具栏的图标都会有提示的，像运行按钮的提示就是Run Figure.我们会看到如下的界面： 那说明，我们保存的.fig 文件的目录不是当前目录，但是没关系啊，我们只要点击“Change Directory ”来改变当前目录。当然，如果你想把当前目录添加到MATLAB 路径也可以，那就点击“Add to Path ”就OK 了。我在这里推荐点击“Change Directory ”，因为没有什么太大必要把其添加到MATLAB 路径中，一般是工具箱需要添加或者我们的函数或程序写完了，而在MATLAB 的命令窗口找不到我们的函数的时候，我们可以将函数或程序所在的目录添加到MATLAB 路径。 总之吧，点那个按钮，要看个人的爱好了。不管点击两个按钮的那一个按钮，都会正确的运行程序的。 我们的程序运行时的样子，是这样的：
程序代码：(代码标记[code]...[/code] ) 1-32是：图形应用篇 33-66是：界面设计篇 67-84是：图形处理篇 85-100是：数值分析篇 实例1：三角函数曲线（1） function shili01 h0=figure('toolbar','none',... 'position',[198 56 350 300],... 'name','实例01'); h1=axes('parent',h0,... 'visible','off'); x=-pi:0.05:pi; y=sin(x); plot(x,y); xlabel('自变量X'); ylabel('函数值Y'); title('SIN( )函数曲线'); grid on 实例2：三角函数曲线（2） function shili02 h0=figure('toolbar','none',... 'position',[200 150 450 350],... 'name','实例02'); x=-pi:0.05:pi; y=sin(x)+cos(x); plot(x,y,'-*r','linewidth',1); grid on xlabel('自变量X'); ylabel('函数值Y'); title('三角函数');
实例3：图形的叠加 function shili03 h0=figure('toolbar','none',... 'position',[200 150 450 350],... 'name','实例03'); x=-pi:0.05:pi; y1=sin(x); y2=cos(x); plot(x,y1,... '-*r',... x,y2,... '--og'); grid on xlabel('自变量X'); ylabel('函数值Y'); title('三角函数'); 实例4：双y轴图形的绘制 function shili04 h0=figure('toolbar','none',... 'position',[200 150 450 250],... 'name','实例04'); x=0:900;a=1000;b=0.005; y1=2*x; y2=cos(b*x); [haxes,hline1,hline2]=plotyy(x,y1,x,y2,'semilogy','plot'); axes(haxes(1)) ylabel('semilog plot'); axes(haxes(2)) ylabel('linear plot'); 实例5：单个轴窗口显示多个图形 function shili05 h0=figure('toolbar','none',... 'position',[200 150 450 250],... 'name','实例05'); t=0:pi/10:2*pi;
实验一特殊函数与图形 一、问题背景与实验目的 二、相关函数（命令）及简介 三、实验内容 四、自己动手 一、问题背景与实验目的 著名的Riemann函数大家都很熟悉了，但是关于它的图像你是否清楚呢？除了最上面那几点，其他都很难画吧？你想不想看看下面那些“挤在一起”的点是怎样分布的呢？还有几何中的马鞍面、单叶双曲面等是怎样由直线生成的，是不是也想目睹一下呢？这些，都离不开绘图． 实际上绘图一直是数学中的一种重要手段，借助图形，往往可以化繁为简，使抽象的对象得到明白直观的体现．比如函数的基本性质，一个图形常可以使之一目了然，非常有效．它虽不能代替严格的分析与证明，但在问题的研究过程中，可以帮助研究人员节约相当一部分精力．此外，它还可以使计算、证明、建模等的结果得到更明白易懂的表现，有时，这比科学论证更有说服力．同时，数学的教学与学习过程也离不开绘图．借助直观的图形，常可以使初学者更容易接受新知识．如数学分析中有不少函数，其解析式着实让人望而生畏，即使对其性质作了详尽的分析，还是感到难明就里；但如果能看到它的图形，再配合理论分析，则问题可以迎刃而解．又如在几何的学习中，会遇到大量的曲线与曲面，也离不开图形的配合． 传统的手工作图，往往费力耗时，效果也不尽理想．计算机恰恰弥补了这个不足，使你可以方便地指定各种视角、比例、明暗，从各个角度进行观察．本实验通过对函数的图形表示和几个曲面（线）图形的介绍，一方面展示它们的特点，另一方面，也将就Matlab软件的作图功能作一个简单介绍．大家将会看到，Matlab 的作图功能非常强大． 二、相关函数（命令）及简介 1．平面作图函数：plot，其基本调用形式： plot(x,y,s) 以x作为横坐标，y作为纵坐标．s是图形显示属性的设置选项．例如：x=-pi:pi/10:pi; y=sin(x); plot(x,y,'--rh','linewidth',2,'markeredgecolor','b','markerfacecolor','g')
MATLAB大作业 一、下料问题的优化设计 1、设计目标 边角总长最短 2、已知参数 下料长度与下料方案 3、MATLAB优化程序及优化结果 （1）取设计变量的初值为x0=[2 6]。首先编写目标函数的m文件Objfun.m，返回x处的函数值f。 （2）运行run.m函数，给初定值，并调用优化函数。 （3）运行结果如下： 即当参数为上述值时，这时所剩余的边角料总长度最短，其值为425cm。 二、键结构优化设计
将键的尺寸b与h分别定为设计变量x1和x2，则X=[x1,x2]T，若按键的体积最小为目标建立目标函数，则目标函数为v=lbh。 1、MATLAB优化程序及优化结果 （1）编写目标函数的m文件myfun.m。 （2）编写非线性约束函数的m文件mycon.m。 （3）在MATLAB命令窗口调用优化程序。 经运行，得到结果如下： （4）优化结果分析。 在强度允许的范围之内，取b=2.5926mm,h=7.7778mm即可满足要求且键的体积最小。 三、单缸内燃机气缸直径和行程的最优化设计
（1）首先编写目标函数的ObjFunct.m文件。 （2）编写返回约束值的函数NonLinConstr.m。 （3）MATLAB主函数代码run.m。 （4）程序运行后，计算结果如下：
四、个人设计的题目： 对长，宽分别为4,6米的长方形铁板，在四个角剪去相等的正方形以制成长方体无盖箱体，问如何剪法使箱体的容积最大？ 解：设剪去正方形的边长为x，则箱体的容积为(6-2x)(4-2x)x 建立无约束优化模型为：miny=(6-2x) (4-2x)x, (0 本科生实验报告 实验课程《MATLAB仿真技术》 学院名称核技术与自动化工程学院 专业名称核工程与核技术 学生姓名杨洲 学生学号201506020219 指导教师杨剑波 实验地点6C802 实验成绩 二〇一六年十月二〇一六年十一月 1.对于典型二阶系统 编写程序并绘制出无阻尼自然振荡ωn=8，阻尼比ξ分别为0.2,0.4,…,1.0,2.0时系统的单位阶跃响应曲线 程序: %典型二阶系统的单位阶跃响应曲线 clc,clear; Wn=8; %无阻尼自然振荡 t=linspace(0,5); for a=0.2:0.2:0.8 %（0<阻尼比a<1,欠阻尼） Wd=Wn.*sqrt(1-a.^2); c=1-exp(-a.*Wn.*t).*(cos(Wd.*t)+a/sqrt(1-a.^2).*sin(Wd.*t)); plot(t,c); hold on end hold on a=1.0; %（阻尼比a=1,临界阻尼） c=1-(1+Wn.*t).*exp(-Wn.*t); plot(t,c); hold on a=2.0; %(阻尼比a=2.0>1,过阻尼） s1=-a.*Wn+sqrt(a.^2-1).*Wn; s2=-a.*Wn-sqrt(a.^2-1).*Wn; c=1-1/(2.*sqrt(a.^2-1)).*(exp(s1.*t)/(a-sqrt(a.^2-1))-exp(s2.*t)/ (a+sqrt(a.^2-1))); plot(t,c); xlabel('t'),ylabel('c(t)'); title('单位阶跃响应曲线') 绘制图形如下： 数字信号处理MATLAB编程作业 姓名：白焱学号：2012001020006 M2.2:代码：（以其中一组实例，其他组相同，其他类似情况不在说明） n1=-10:10; xn=5*cos(1.5*pi*n+0.75*pi)+4*cos(0.6*pi*n)-sin(0.5*pi*n); stem(xn); xlabel('n'); ylabel('xe[n]'); 结果： M2.9代码： clear;clc; x = input('Type in the reference sequence = '); y = input('Type in the second sequence = '); % Compute the correlation sequence n1 = length(y)-1; n2 = length(x)-1; r = conv(x,fliplr(y)); k = (-n1):n2'; stem(k,r); xlabel('Lag index'); ylabel('Amplitude'); v = axis; axis([-n1 n2 v(3:end)]); 结果： X[n]自相关序列 y[n]自相关序列 w[n]自相关序列 X[n]与y[n]互相关序列 X[n]与w[n]互相关序列 M3.1 代码： rr=input('value of r ='); xita=input('value of xita ='); fem(1)=0; fem(2)=-2*rr*cos(xita); fem(3)=rr*rr; num = [1]; den = fem; % Compute the frequency response w = 0:pi/(k-1):pi; h = freqz(num, den, w); % Plot the frequency response subplot(2,2,1) plot(w/pi,real(h));grid title('Real part') xlabel('\omega/\pi'); ylabel('Amplitude') subplot(2,2,2) plot(w/pi,imag(h));grid title('Imaginary part') xlabel('\omega/\pi'); ylabel('Amplitude') subplot(2,2,3) plot(w/pi,abs(h));grid title('Magnitude Spectrum') xlabel('\omega/\pi'); ylabel('Magnitude') subplot(2,2,4) plot(w/pi,angle(h));grid title('Phase Spectrum') xlabel('\omega/\pi'); ylabel('Phase, radians') 现代设计方法-工程优化理论、方法与设计 姓名 学号 班级 研 问题 ： 某厂向用户提供发动机，合同规定，第一、二、三季度末分别交货40台、60台、80台。每季度的生产费用为 （元），其中x 是该季生产的台数。若交货后有剩余，可用于下季度交货，但需支付存储费，每台每季度c 元。已知工厂每季度最大生产能力为100台，第一季度开始时无存货，设a=50、b=0.2、c=4，问工厂应如何安排生产计划，才能既满足合同又使总费用最低。讨论a 、b 、c 变化对计划的影响，并作出合理的解释。 问题的分析和假设： 问题分析：本题是一个有约束条件的二次规划问题。决策变量是工厂每季度生产的台数，目标函数是总费用（包括生产费用和存储费）。约束条件是生产合同，生产能力的限制。在这些条件下需要如何安排生产计划，才能既满足合同又使总费用最低。 问题假设： 1、工厂最大生产能力不会发生变化； 2、合同不会发生变更； 3、第一季度开始时工厂无存货； 4、生产总量达到180台时，不在进行生产； 5、工厂生产处的发动机质量有保证，不考虑退货等因素； 6、不考虑产品运输费用是否有厂家承担等和生产无关的因素。 符号规定： x1——第一季度生产的台数； x2——第二季度生产的台数； 180-x1-x2——第三季度生产的台数； y1——第一季度总费用； y2——第二季度总费用； y3——第三季度总费用； y ——总费用（包括生产费用和存储费）。 ()2bx ax x f += 建模： 1、第一、二、三季度末分别交货40台、60台、80台； 2、每季度的生产费用为 （元）； 3、每季度生产数量满足40 ≤x1≤100，0≤x2≤100，100≤x1+x2 ≤180； 4、要求总费用最低，这是一个目标规划模型。 目标函数: y1 2111x b x a Z ?+?= y2()4012222-?+?+?=x c x b x a Z y3()()()100180180212 21213-+?+--?+--?=x x c x x b x x a Z y x x x x x x Z Z Z Z 68644.04.04.0149201212 221321--+++=++= 40≤x1≤100 0≤x2≤100 100≤x1+x2≤180 ()2bx ax x f +=matlab简单实例作业
数字信号处理MATLAB编程作业
实例matlab 非线性规划 作业
相关主题
文本预览