mathematica讲座 ----徐安农教授 .txt

Matlab在化学工程中的应用

b={-850,-700,-100,-900};  
LinearProgramming[c,A,b]  
(* 此命令只给出决策变量。*)  
  
线性约束条件下的非线性规划问题  
线性逼近法(FW法)  
模型: (NLP)minf(x)  
S.t.: E={x|AX>=b,X>=0}  
执行算例4 求解非线性规划问题  
用线性逼近法求解非线性规划问题  
目标函数 f(x1,x2)=(x1-1)^2+(x2-2)^2  
约束条件 0<=x1<=2,0<=x2<=3  
下面是第一次迭代  
Clear[a,b,c,d,s,e,pf]  
f[x1_,x2_]:=(x1-1)^2+(x2-2)^2;  
gradf={D[f[x1,x2],x1],D[f[x1,x2],x2]};  
c={0.7,1.25}; (*C即基可行解X0*)  
s=gradf/.{x1->Part[c,1],x2->Part[c,2]} (*求X0点梯度*)  
{-0.6, -1.5}  
x1=.; (*清除X1,X2的值*)  
x2=.;  
p[u_,v_]:=s.{u,v}  
a=ConstrainedMin[p[x1,x2],{x1<=2,x2<=3},{x1,x2}]  
(*求出最优解Y0*)  
b={x1,x2}/.Part[a,2];  
e=b-c;  
pf=s.e;  
If[Abs[pf]>0.01,  
g=c+d*e;  
t[d_]:=f[Part[g,1],Part[g,2]];  
w=FindMinimum[t[d],{d,1,0,1}];  
c=c+(d/.Part[w,2])e;  
]  
Print["c1=",c]  
(* 得到初始点c1,将其替换c,运算后的结果继续代替C,直到w的绝对  
值小于0.01为止。*) *)  
  
  
MATHEMATICA讲座第九讲  
表处理函数  
　Mathematica的表处理功能非常强大，这成为Mathematica系统的一个显著特点，  
在符号处理、矩阵运算等方面具有重要作用,在一般的程序设计中亦是必不可少的。 
1、Table[表达式，循环描述]：按循环描述生成具有表达式所描述的性质的表；  
算例：  
Table[i2+1，{i，1，5}]  
2、(*NestList[函数f，初值x，递推次数n]：求值生成一个由n+1元素的表：  
　　　　　　　　　　　　　{f（x），f（f（x）），…，f…f（x）}　*)  
算例：  
h[x_]:=x^2;  
NestList[h,2,4]  
3、Part[表，n]（或表[[n]]）：取出表的第n个元素；  
算例：  
Clear[a,b,c,d,e];  
ww={a,b,c,d,e};  
ww[[3]]  
Part[ww,2]  
4、Take[表，整数n]：取出表的前n个元素做成一个表，  
Drop[表，n]的作用正好与Take相反，是取出后表的前n个元素剩余的元素作成的一个  
表；  
算例：  
Clear[a,b,c];  
Take[{1,3,5,a,b,c},4]  
Drop[{1,3,5,a,b,c},2]  
5、Count[表，表达式]：求表达式在表的第一层出现的次数；  
算例：  
Clear[a,b,c,d];  
Count[{a,b,c,a,d,a,{a,c}},a]  
6、Position[表，表达式]；找出表达式在表中出现的位置；  
算例：  
Position[{1，b，a，{a，b}}，a]]  
7、Prepend[表，表达式]：把表达式放在原表所有元素的前面；  
Append[表，表达式]：表达式放在原表最后  
算例：  
Prepend[{a,b,c,s},f]  
Append[{a，b，c，d}，e]  
8、Insert[表，表达式，整数n]：表达式插入原表第n个位置  
算例：  
Insert[{1，2，3，4}，a，3]  
9、Join[表1，表2，…]：得到一个由这n个表的元素顺序连接起来构成的表；  
Union[表1，表2，…]：与Join类似，只是在作为结果的表里删除了重复元素，并且  
将表的元素按一种内定的次序重新排序；  
算例：  
Join[{a,b,c},{a,b}]  
Union[{a，e，f}，{a，b，e，f}]  
10、Reverse[表]：求出原表反序的表；  
算例：  
Reverse[{1,2,3,4}]  
11、Transpose[表]：这里的表应多于一层，求出原表第一层和第二层元素交换得的表，如果表为矩阵即得矩阵的转置；  
算例：  
Transpose[{{1,2,3},{4,5,6}}]  
12、 Deletecases[表，表达式]：删除表中与表达式相同的元素；  
算例：  
Deletecases[{10,8,9,7},10]  
13、Ｆlatten[表，n]：求出表上面n层抹平后得到的表；  
算例：  
Flatten[ { { a，b，{c，d，{ e，f } } }，g }，2 ]  
14、 Det[矩阵]：求矩阵的行列式；  
Inverse[矩阵]：求矩阵的逆；  
算例：  
A={{1,2,3},{4,5,6},{7,8,8}};  
Det[A]  
Inverse[A]  
15、Intersection[表1，表2，…]：将这些表作为集合求交集；  
算例：  
Clear[a,b,c,d,e,g];  
Intersection[{a,b,c,d},{a,b,e,g},{b,c,d}]  
16、Complement[表1，表2，表3，…]：求表1，表2，表3，…相对于 表1的补集；  
算例：  
Clear[a,b,c,d,e];  
Complement[{a，b，c，d}，{a，c}，{a，e}]  
17 Map[函数，表达式]：将函数作用到表达式的第一层的每一个元素上，得到由这样作用的结果构成的表达式；  
算例：  
Clear[a,b,c];  
Map[#Λ2& ，{a，b，c}]  
18、Apply[函数，表]：把表作为函数的参数表求值得到的表达式  
算例：  
Clear[a,b,c]  
ff[x_,y_,z_]:=x+y^2+z^3;  
Apply[ff,{a,b,c}]  
  
  
MATHEMATICA第十讲 
程序设计初步  
在这一讲里我们将介绍 Mathematica的程序设计的初步知识。  
系统自身定义了几百个函数，这是MATHEMATICA之所以功能强大的缘故，很好的利用这些函数可以完成许多方面的工作，如 Fit Plot等。这里学习一个很有用的指令“？”。  
  
在一个函数前面打 ？可以得到有关这个函数的说明。  
例如：  
?Plot  
在一个字母前加？，系统给出以这个字母开头的所有函数的列表。  
?P*  
??Plot 给出关于Plot的进一步的信息  
函数的复合  
利用系统内部函数可以构造出更多更复杂的函数。 
例 先求出SIN(X)的台劳展开式，再截取它的前九项。  
Series[Sin[x],{x,0,9}]  
Normal[Series[Sin[x],{x,0,9}]]  
  
二 自定义函数  
f[x_]:=表达式  
如果输入的自变量是表,则用:  
f[x_List]:=表达式  
一般地f[x_patten]:=表达式,patten表示自变量的模式,如:  
f[x_Integer]:=表达式,自变量应为整数.  
同理f[x_,y_]:=表达式,定义二元函数.  
例  
f[x_]:=x^3-x+1  
Plot[f[x],{x,-2,2}]  
Print["f(0)=",f[0]]  
例 定义一个画半径为r的圆的函数myPlot[r_]  
myPlot[r_]:=ParametricPlot[{r*Cos[t],r*Sin[t]},{t,0,2*Pi},  
AspectRatio->Automatic]  
myPlot[1]  
Map[myPlot,{1,2,3,4,5}]  
Show[%]  
  
三 循环结构  
1 While[条件，表达式]  
称为当循环结构。其中“条件”是一个逻辑表达式，先对条件求值，如得到True,  
则求值它的表达式部份，然后重复上述过程，直到条件不满足为止，循环结束。  
例 求平方小于100的最大的整数  
clear[x];  
x=0;  
While[x^2<=100,x=x+1];  
Print["x=",x-1]  
例 用黄金分割法求方程的根  
eps=10^(-6);a=-1.;b=0.;  
f[x_]:=x^2+2x+2;  
x1=b-0.618(b-a);x2=a+0.618(b-a);  
f1=f[x1];f2=f[x2];  
While[x2-x1>eps,  
If[f2>f1,b=x2;x2=x1;x1=b-0.618(b-a);f2=f1;f1=f[x1],  
a=x1;x1=x1;x2=a+0.618(b-a);f1=f2;f2=f[x2]  
]  
];  
Print["最优解x*=",(x1+x2)/2," 最优值f(x*)=",f[(x1+x2)/2]]  
2.For[初始表达式，条件，步进表达式，表达式]  
步进表达式用于对循环控制变量作步进赋值  
n++表示将n的值增加一个单位;  
n+=2表示将n的值增加两个单位;  
例 求出前10个自然数的和L与前10个自然数的乘积S  
L=0;S=1;  
For[i=1,i<=10,i++,  
L=L+i;  
S=S*i;  
]  
Print["L=",L]  
Print["S=",S]  
例 计算1+1/3+1/5+...+1/10  
S=0;  
For[i=1,i<=10,i+=2,  
S=S+1/i;  
]  
Print["S=",N[S,10]]  
  
例 求Fibonacci数列的前40项  
fibonacci={1,1};  
For[i=1,i<=38,i++,  
fibonacci=Append[fibonacci,fibonacci[[i]]+  
fibonacci[[i+1]]]  
]  
Print["fibonacci[[40]]=",fibonacci[[40]]]  
例 求自然对数的底e的近似值  
e=1;t=1;  
For[k=1,k<=10,k++,  
t=k*t;  
e=e+1/t  
];  
N[e,10]  
3 Do[表达式，循环描述]  
例 将表达式t=t^2+1从t=1执行3次  
Clear[t]  
t=1;  
Do[t=t^2+1,{3}]  
Print["t=",t]  
t=26  
例 计算翻倍  
Clear[t]  
t=5;  
Do[t=2*t,{n,2,5}]  
Print["t=",t]  
例  
Clear[t]  
Do[Print[t^2],{t,4}]  
1  
例：牛顿迭代法求方程x-cos(x)=0的实根  
Plot[x-Cos[x],{x,-10,10}]  
t=3;  
Do[Print[t];t=N[t-(t-Cos[t])/(1+Sin[t])],{20}  
]  
4 FixedPoint[函数，表达式]  
将已给函数一直复合到不动点，如果存在的话。  
  
例 求证线性函数0.2*x+5无论从哪一个点出发经过多次迭代都将收敛于同一个数。  
f[x_]:=0.2 x+5;  
FixedPoint[f,0]  
5 表达式//.替换规则  
例 用t+1代换x^2+2*x-1之中的x直到不变为止。  
Clear[t,x]  
x^2+2*x-1//.x->t+1  
6 Nest[函数，表达式，整数n] 将表达式代入函数作用给定的次数。  
例  
f[x_]:=x^2;  
Nest[f,x^2+x+1,3]  
  
四 分枝结构  
1 If [条件，表达式]  
当条件为True时，求表达式的值，当条件为False时，返回Null.  
If[条件，表达式1，表达式2]  
当条件为True时，求表达式1的值，当条件为False时，求表达式2的值。  
If[条件，表达式1，表达式2,表达式3]  
当条件为True时，求表达式1的值，当条件为False时，求表达式2的值,  
当条件表达式求不出True或False时,以表达式3的值作为结果  
例  
Clear[x,y];  
x=10  
If[x!=Integer,y=1,y=2,y=3];  
y  
例  
Clear[x,y];  
x=0.4;  
If[x!=Integer,y=1,y=2,y=3];  
y  
例  
Clear[x,y,tt];  
x=tt;  
If[x!=Integer,y=1,y=2,y=3];  
y  
例 观察数列的前50项  
求极限 xn=Lim(1+1/2^2+1/3^2+...+1/n^2)  
xn1=0;m=4;xn=1.;  
For[i=2,i<=50,i++,  
If[Abs[xn-xn1]>10^(-m),  
xn1=N[xn,10];  
xn=xn+1/i^2;  
Print[i," ",xn1," ",xn," "]  
]  
]  
例 求100以内的素数  
w={};  
For[i=1,i<=100,i++,  
If[PrimeQ[i],w=Append[w,i]]  
];  
w  
Length[w] (*给出100内的素数个数*)  
例 构造一个函数求x以内的所有素数  
Prn[x_]:=Module[{i,w={}},  
For[i=1,i<=x,i++,  
If[PrimeQ[i],AppendTo[w,i]]];w  
]  
Timing[Prn[1000]]  
(*其中Timing[]是显示计算的时间。*)  
2 Which[条件1，表达式1，条件2，表达式2,...]  
例 用Which定义分段函数  
f[x_]:=Which[x>0,x+1,x<0,exp[x]];  
Print["f[2]=",f[2]]  
Print["f[-1]=",f[-1]]  
Print["f[0]=",f[0]]  
例 按分数段统计一个教学班某次考试成绩学生人数  
成绩数据  
76，58，84，32，91，95，94，88，78，83，82，67，63，  
69，74，77，100，78，81，92，95，67，49，90，53，67，  
67，87，78，94  
data={76,58,84,32,91,95,94,88,78,83,82,67,63,69,74,  
77,100,78,81,92,95,67,49,90,53,67,67,87,78,94};  
Print["不及格人数为",Length[Select[data,#<60&]]]  
Print["60~69分的人数为",Length[Select[data,#>=60&&#<70&]]]  
Print["70~79分的人数为",Length[Select[data,#>=70&&#<80&]]]  
Print["80~89分的人数为",Length[Select[data,#>=80&&#<90&]]]  
Print["90分以上的人数为",Length[Select[data,#>90&]]]  
3 Switch[判别表达式，模式1，表达式1，模式2，表达式2,...]  
例 构造一个函数，当x被3整除时，函数值为a，当x关于3的模为1时，函数值为b,当x  
关于3的模为2时，函数值为c.  
h[x_]:=Switch[Mod[x,3],0,a,1,b,2,c]  
h[4]  
h[5]  
h[6]  
  
五 转向结构  
Break[] 用于For、While（Do不能使用）语句中,表示结束包含这个break表达式  
的最近循环，且以Null作为该结构的值。  
Continue[] 立即结束包含这个Continue的循环,执行下一次循环。  
Return 从函数的求值过程中退出(例句见十四讲)  
Return[] 退出函数的求值,以Null作为当前函数值  
Return[expr] 退出函数的求值,以expr作为当前函数值  
典型的使用结构为:  
While(或For)[...,  
...,...;  
If[...,Continue[]];  
...,...;  
If[...,...,Break[]];  
...,...  
]  
  
例： 求20至300之间不能被5整除的所有数  
w={};  
For[i=50,i<=300,i++,  
If[Mod[i,5]==0,Continue[]];w=Append[w,i];  
];  
w  
例 求20至500之间不能被8整除的第十个数  
w={};  
For[i=50,i<=300,i++,  
If[Mod[i,8]!=0,w=Append[w,i]];If[Length[w]==10,Break[]]  
];  
i  
Goto[标志]和Lable[标志]用于在复合表达式中实现执行的控制转移.  
Lable[标志]是一个位置的标识，本身什么也不做。  
Goto[标志]将执行立即转移到具有同一标志的Lable位置去。  
典型使用方式  
Module[...,  
......;  
Lable[a];  
......;  
If[...,Goto[a]];  
......]  
例 Goto的作用在于越过某些语句去执行指定的语句  
x=1;y=3;  
If[y>x,Goto[aa]];  
y=78;  
Label[aa];  
y