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

Matlab在化学工程中的应用

Simplity[f[x]] 将表达式化简  
执行算例  
Integrate[1/Sqrt[1-x^2],x]  
执行算例  
Integrate[1/Sqrt[1-x^2],{x,-1,1}]  
执行算例  
Integrate[(x+1)/(x^2+3x+5),x]  
D[%,x](*验证之*)  
Simplify[%]  
执行算例  
NIntegrate[1/Sqrt[1-x^2],{x,-1,1}]  
Plot[Sin[x]/x,{x,-2Pi,2Pi}]  
NIntegrate[Sin[x]/x,{x,0.000001,Pi}]  
  
四.幂级数展开  
Series[f[x],{x,x0,n}] 将函数f(x)在x0点展开n阶  
Normal[%] 去掉幂级数的余项尔后可以求值或画图形  
执行算例  
S=Series[Sin[x],{x,0,10}]  
S1=Normal[S]  
a1=Plot[Sin[x],{x,0,2Pi},PlotStyle->RGBColor[1,0,0] ]  
a2=Plot[S1,{x,0,2Pi}, PlotStyle->RGBColor[0,0,1]]  
Show[a1,a2,AspectRatio->Automatic]  
  
  
MATHEMATICA讲座第六讲  
线性方程组的表达方式和解法  
  
一.向量和矩阵的输入  
  
1.Range[n]  
Range[n,m]  
Range[n,m,h]  
执行算例  
Range[5]  
Range[1,10,2] 
2.Table[(1/2)^n,{n,0,10}] (*由通项构造表*)  
Table[{f1(n),F2(n)},{n,n1,n2,h}]  
执行算例  
Table[(1/2)^n,{n,0,10}]  
A=Table[Normal[Series[Sin[x],{x,0,n}]],{n,1,11,2}]  
(*Sin[x]的六个正规化后的幂级数展开式的表*)  
执行算例  
Table[{x,x^2},{x,-1.0,1.0,0.2}]  
Table[Expand[(x^i+i*x)^2],{i,2,5}]  
Table[Mod[n,2],{n,0,17}]  
Table[Sum[x^i,{i,0,n}],{n,1,5}]  
Table[Random[],{10}]  
Table[Random[Real,{1,10}],{10}]  
3.Array[函数,n] (*由函数表达式构造表*)  
Array[函数，{n1,n2,n3,...}]  
执行算例  
Array[Exp,5]（*等价于Table[Exp[x],{x,5}]*）  
Array[Mod,{10,10}]  
（*等价于Table[Mod[n,m]{n,0,10},{m,0,10}]*）  
4.NestList[f[n],n0,k] 用递推公式建立表元素  
Clear[n，n0,rnt,fnt,t1,t2]  
t1=(Sqrt[5]+1)/2  
t2=(1-Sqrt[5])/2  
fnt=Table[(t1^(n+1)-t2^(n+1))/Sqrt[5],{n,0,40}]//N  
rnt=Table[fnt[[n-1]]/fnt[[n]],{n,2,12}]  
5.向量与矩阵的标准输入法  
A={x1,x2,x3,...}  
A={{a11,a12,a13},{a21,a22,a23},{a31,a32,a33}}  
ColumnForm[{a1,a2,...,an}]把向量用列方式输出  
MatrixForm[A] 用矩阵方式显示  
IdentityMatrix[n] 生成n阶单位阵  
DiagonalMatrix[{a11,a22,...,ann}]生成对角阵  
执行算例  
A1={1,2,3}  
ColumnForm[A1]  
A2=DiagonalMatrix[{1,2,3,4,5}]  
MatrixForm[A2]  
IdentityMatrix[3]  
DiagonalMatrix[{1,2,3,4,5}]  
MatrixForm[%]  
  
二.行列式  
Det[A]  
执行算例  
A={{1,2,3},{4,5,6},{7,8,9}}  
MatrixForm[A]  
Det[A]  
  
三.矩阵求逆,求特征值，特征向量  
Inverse[A]  
Eigenvalue[A]  
Eigenvector[A]  
执行算例  
Clear[A]  
A={{4,6,0},{-3,-5,0},{-3,-6,1}}  
Eigenvalues[A]  
Eigenvectors[A]  
  
四.恰定方程求解  
问题1 x1+6x2+36x3=104  
x1+10x2+100x3=160  
x1+20x2+400x3=370  
程序  
Clear[A1,b]  
A1={{1,6,36},{1,10,100},{1,20,400}};  
b={104,160,370};  
LinearSolve[A1,b](*求方程组的解*)  
X=Inverse[A1].b (*用求逆矩阵方法求解*)  
  
五.欠定方程求解  
问题2 2x1+x2-x3+x4=1  
x1+2x2+x3-x4=2  
x1+x2+2x3+x4=3  
程序  
A2={{2,1,-1,1},{1,2,1,-1},{1,1,2,1}}  
b2={1,2,3};  
X={x1,x2,x3,x4};  
Solve[A2.X==b2,{x1,x2,x3,x4}]  
  
  
MATHEMATICA讲座第七讲  
函数的插值  
  
一.拉格朗日插值  
L={List}  
InterpolatingPolynomial[L,x] 
执行算例1 两点线性插值  
L={{0,0.3},{0.2,0.45}}  
I=InterpolatingPolynomial[L,x]  
执行算例2 三点抛物插值  
L1={{0,0.3},{0.2,0.45},{0.4,0.15}}  
I1=InterpolatingPolynomial[L1,x]  
执行算例3 多点拉格朗日插值  
L2={{0,0.3},{0.2,0.45},{0.3,0.47},  
{0.52,0.50},{0.64,0.38},{0.7,0.33},{1.0,0.24}}  
I2=InterpolatingPolynomial[L2,x]  
Plot[%,{x,-0.25,1.05}]  
执行算例4 作正弦在0，P上五点插值函数图形  
g0=Plot[Sin[x],{x,0,Pi}]  
L=Line[Table[{x,Sin[x]},{x,0,Pi,Pi/4}]]  
g=Graphics[L]  
Show[g0,g]  
sinAp[n_]:=Graphics[{Line[Table[{x,Sin[x]},  
{x,0,Pi,Pi/(n+1)}]]}]  
sinAp[2]  
Show[g0,%]  
  
二.龙格现像演示  
L=Table[{x,1/(1+25*x^2)},{x,-1,1,0.2}]  
a=InterpolatingPolynomial[L,x]  
b=Plot[1/(1+25*x^2),{x,-1,1},  
PlotStyle->{RGBColor[1,0,0]}]  
c=Plot[a,{x,-1,1}]  
Show[b,c]  
  
三. 两点三次Hermite插值  
执行算例5  
Clear[x,y,h0,h1,H0,H1]  
x1={0,1};y1={1,2};m={1/2,1/2};  
h0[x_]:=(1+2*x)*(x-1)^2;  
h1[x_]=(1-2(x-1))*(x/(x-1))^2  
H0[x_]=x*(x-1)^2  
H1[x_]=(x-1)*x^2  
H[x_]=y1[[1]]*h0[x]+y1[[2]]*h1[x]  
+m[[1]]*H0[x]+m[[2]]*H1[x]  
%/.{x->0.55}  
  
四. N+1个节点的2N+1次Hermite插值  
执行算例6  
Clear[x0,y,bb,w,w1,w2,L,h,H,Hm]  
x0={0.4,0.5,0.6,0.7,0.8}  
y=Table[Log[x],{x,0.40,0.80,0.10}]  
m=Table[1/x,{x,0.40,0.80,0.10}]  
bb[x_]=InterpolatingPolynomial[y,x]  
Simplify[bb[x]]  
bb[0.55]  
w[x_]=(x-x0[[1]])*(x-x0[[2]])*(x-x0[[3]])*(x-x0[[4]])*(x-x0[[5]])  
w1[x_]=D[w[x],x];  
Simplify[w1[x]]  
w2[x_]=D[w[x],{x,2}];  
Simplify[w2[x]]  
For[ i=1,i<=5,i++,  
L[i_,x_]:=w[x]/((x-x0[[i]])*w1[x0[[i]]])];  
h[i_,x_]:=(1-w2[x0[[i]]]*(x-x0[[i]])/w1[x0[[i]]])*L[i,x]^2;  
H[i_,x_]:=L[i,x]^2*(x-x0[[i]]);  
]  
Hm[x_]=Sum[y[[i]]*h[i,x]+m[[i]]*H[i,x],{i,1,5,1}];  
Simplify[Hm[x]]  
Hm[0.55]  
  
拟合  
  
一.一元线性拟合  
执行算例  
b2={{100,45},{110,51},{120,54},{130,61},{140,66},  
{150,70},{160,74},{170,78},{180,85},{190,89}}  
fp=ListPlot[b2,PlotStyle->{PointSize[0.03],  
RGBColor[0,0,1]}]  
ft1=Fit[b2,{1,x},x]  
gp=Plot[ft1,{x,100,190},PlotStyle->{RGBColor[1,0,0]}];  
Show[fp,gp]  
  
二.抛物线拟合  
执行算例  
B=Table[Prime[n],{n,20}]  
t1=ListPlot[B,PlotStyle->{RGBColor[0,1,1],  
PointSize[0.04]}]  
f=Fit[B,{1,x,x^2},x]  
t2=Plot[f,{x,0,20}]  
Show[t1,t2]  
  
三.多项式拟合  
执行算例  
data={{0,1.2},{1,1.4},{2,1.3},{3,1.5},{4,1.3},{5,1.3},{6,1.1}};  
t2=ListPlot[data,PlotStyle->{PointSize[0.05],  
RGBColor[1,0,0]}]  
fx=Fit[data,{1,x,x^2,x^3,x^4,x^5},x]  
t1=Plot[fx,{x,0,6},PlotStyle->RGBColor[1,0,1]]  
Show[t1,t2]  
执行算例  
b3={{1,4},{2,6.4},{3,8.0},{4,8.4},{5,9.28},{6,9.5},  
{7,9.7},{8,9.86},{9,10.0},{10,10.2},{11,10.32},  
{12,10.42},{13,10.5},{14,10.55},{15,10.58},  
{16,10.6}}  
gp=ListPlot[b3,PlotStyle->{RGBColor[0,1,0],  
PointSize[0.04]}]  
ft2=Fit[b3,Table[x^i,{i,0,4}],x]  
fp=Plot[ft2,{x,0,17},PlotStyle->{RGBColor[1,0,0]}]  
Show[gp,fp]  
  
四.非线性拟合---指数拟合  
执行算例 求一个经验函数,型如  
x 1 2 3 4 5 6 7 8  
y 15.3 20.5 27.4 36.5 49.1 65.5 87.8 117.6  
程序  
b4={{1,15.3},{2,20.5},{3,27.4},{4,36.6},{5,49.1},  
{6,65.5},{7,87.8},{8,117.6}};  
gb4=ListPlot[b4,PlotStyle->{RGBColor[0,0,1],  
PointSize[0.05]}]  
y4=Table[b4[[i,2]],{i,1,8}];  
ly4=Log[y4];  
fy4=Fit[ly4,{1,x},x]  
s[x_]:=Exp[fy4]  
ty=Plot[s[x],{x,1,8},Axes->{2,60},  
AspectRatio->1,  
PlotStyle->{RGBColor[1,0,0]},  
PlotRange->{10,120}]  
Show[ty,gb4, Axes->{2,60},AxesLabel->{"x","y"},  
AspectRatio->1,  
Ticks->{{1,2,3,4,5,6,7,8},{0,30,60,90,120}}]  
执行算例 求一个经验函数,型如y=a*exp(-bx)与所给数据拟合  
x 0.4 0.5 0.6 0.7  
y 1.75 1.34 1.00 0.74  
程序  
Clear[fx,fy,biao,nb.ft,ft1,t1]  
fx[x_]:=x  
fy[y_]:=Log[y]  
biao={{0.4,1.75},{0.5,1.34},{0.6,1.00},{0.7,0.74}};  
nb=Table[{fx[biao[[i,1]]],fy[biao[[i,2]]]},{i,1,4}];  
(*拟合方程*)  
ft=Fit[nb,{1,x},x]  
ft1=Exp[ft]  
(*拟合曲线*)  
t1=Plot[ft1,{x,0,1.0},AxesLabel->{"x","y"},  
PlotStyle->{RGBColor[1,0,0]}]  
t2=ListPlot[biao,PlotStyle->{RGBColor[0,0,1],PointSize[0.04]}]  
Show[t1,t2,PlotRange->{0,2}]  
  
五.用正交多项式作拟合 
[0,1]区间上的勒让得多项式 
(*定义勒让得函数(n=10)*)  
Clear[x,t,s]  
n=10  
P0[x_]:=1  
P1[x_]:=1-2x/n  
P2[x_]:=1-6 x/n+6 x*(x-1)/(n*(n-1))  
P3[x_]:=1-12 x/n+30x*(x-1)/(n*(n-1))-20 x*(x-1)*(x-2)/(n*(n-1)*(n-2))  
P4[x_]:=1-2 x+x*(x-1)-140x*(x-1)*(x-2)/720+70x*(x-1)*(x-2)*(x-3)/(10*9*8*7)  
P5[x_]:=1-30x/n+210x*(x-1)/(n*(n-1))-560x*(x-1)*(x-2)/(n*(n-1)*(n-2))  
+630x*(x-1)*(x-2)*(x-3)/(n*(n-1)*(n-2)*(n-3))  
-252 x*(x-1)*(x-2)*(x-3)*(x-4)/(n*(n-1)*(n-2)*(n-3)*(n-4))  
  
  
(*输入初始数据*)  
t={0,5,10,15,20,25,30,35,40,45,50};  
y={0,1.27,2.16,2.86,3.44,3.87,4.15,4.37,4.51,4.60,4.66};  
(*做变量替换*)  
x=t/5;  
(*计算各多项式在节点处的值*)  
A={P0[x],P1[x],P2[x],P3[x],P4[x],P5[x]}  
(*计算每一行元素平方的和*)  
s=Table[0,{i,1,6}];  
For[k=1,k<=6,k++,  
s[[k]]=0;  
For [i=1,i<=11,i++,  
s[[k]]=s[[k]]+A[[k,i]]^2]  
]  
N[s,6]  
(*计算Pk(xi)*yi*)  
r=Table[0,{i,1,6},{j,1,11}]  
For[i=1,i<=11,i++,  
r[[1,i]]=b0[[i]]*y[[i]]  
]  
For[i=1,i<=11,i++,  
r[[1,i]]=b0[[i]]*y[[i]]  
]  
For[i=1,i<=11,i++,  
r[[1,i]]=b0[[i]]*y[[i]]  
]  
For[i=1,i<=11,i++,  
r[[1,i]]=b0[[i]]*y[[i]]  
]  
For[i=1,i<=11,i++,  
r[[2,i]]=b1[[i]]*y[[i]]  
]  
For[i=1,i<=11,i++,  
r[[3,i]]=b2[[i]]*y[[i]]  
]  
For[i=1,i<=11,i++,  
r[[4,i]]=b3[[i]]*y[[i]]  
]  
For[i=1,i<=11,i++,  
r[[5,i]]=b4[[i]]*y[[i]]  
]  
For[i=1,i<=11,i++,  
r[[6,i]]=b5[[i]]*y[[i]]  
]  
r  
  
  
MATHEMATICA讲座第八讲  
线性规划与非线性规划  
  
模型 min(or max) f=c'X=c1x1+c2x2+...+cnxn  
s.t. MX>=b,X>=0  
  
命令格式  
ConstrainedMax[f,{inequalities},{x,y,...}]  
ConstrainedMin[f,{inequalities},{x,y,...}]  
LinearProgramming[c,M,b] 
执行算例1 求解线性规划问题  
max f=2x+3y  
s.t. x+2y<=8  
0<=x<=4,0<=y<=3  
ConstrainedMax[2*x+3*y,{x+2 y<=8,x<=4,y<=3},{x,y}]  
第一个参数是目标函数的最优值，第二个参数是决策变量的取值。  
执行算例2  
min f=x+2y+3z  
s.t. 2x-y =1  
x +z=1  
x,y,z>=0  
(*解法一*)  
ConstrainedMin[x+2y+3z,{2x-y==1,x+z==1},{x,y,z}]  
(*解法二*)  
c={2,-3}  
M={{-1,-1},{1,-1},{1,0}}  
b={-10,2,1}  
q=LinearProgramming[c,M,b]  
执行算例3 求解线性规划问题  
min -f =- 0.40x1-0.28x2-0.32x3-0.72x4-0.64x5-0.61x6  
s.t. -0.01x1-0.01x2-0.01x3-0.03x4-0.03x5-0.03x6>=-850  
-0.02x1 -0.05x4 >=-700  
-0.02x2 -0.05x5 >=-100  
-0.03x3 -0.08x6>=-900  
x1,x2,x3,x4,x5,x6>=0  
解：  
c={-0.4,-0.28,-0.32,-0.72,-0.64,-0.6};  
A={{-0.01,-0.01,-0.01,-0.03,-0.03,-0.03},  
{-0.02,0,0,-0.05,0,0},{0,-0.02,0,0,-0.05,0},  
{0,0,-0.03,0,0,-0.08}};