Mathematica语言.doc

数学模型讲义 林健良附录 Mathematica 软件简介Mathematica是一个功能强大的数学软件.它集数值计算、符号运算,绘图功能于一身,堪称众多数学软件中的佼佼者.加之其语法规则简单,操作使用方便,深受广大科技工作者的喜爱,得到广泛的使用.数学函数和常数 Mathematica提供了大量的数学函数，给运算带来很大方便下面列出一些常用的函数 函 数 形 式功 能Sqrtx平方根Expx指数函数exLogx,Logb,x对数函数Lnx,LogbxSinx,Cosx,Tanx三角函数ArcSinx, ArcCosx, ArcTanx反三角函数Sinhx,Coshx,Tanhx双曲函数n!,n!阶乘,双阶乘Binomialn,m组合数CnmAbsx绝对值Signx符号Roundx四舍五入取整Floorx取不超过x的最大整数Modn,mn/m的余数Random，RandomInteger,m,n，RandomReal,a,b，均匀分布随机数Maxx,y,Minx,y,最大值，最小值Sumai,i,imin,imax,求和Productai, i,imin,imax求积Absz,Argz模，辐角Rez,Imz实部，虚部Conjugatez共轭复数注： Mithematica提供的函数，其名称中的字母大小写是固定的（特别开头字母均为大写），不得误用；函数的自变量以方括号 括起来 Mathemaica还提供了许多数学常数，下面列出了一些常数（均以大写字母开头） Pi -; E-e I-; Infinity- 函数和常数均可参与运算，下面是一些运算的例子 In l： Pi2 Out 12 In2：N Pi,11 Out23.1415626535 In3：LogE8 Out38 In4：SinSqrt1/6 Out4=1/2 用户不仅可以使用Mathemaica提供的函数和常数，还可以自定义函数和常数方法如下： 形式功能 fx_:= expr-定义函数 f fx_,y_:=exp r-定义多变量的函数f ？f-显示函数的定义 Clearf-清除 f的定义 xvalue-给变量x赋值 x清除变量x的值注：定义函数时，在等式左端的方括号中的变量必须跟随下到线符号“”；定义的函数或变量的名称不要使用大写字母开头，以免和Mathemaica的函数或常数混淆例： In1：fx_:=x5；fx_,y_:=Sqrtx2+y2；z3；其中输入语句后的分号“；”表示不显示输出结果，定义了函数、变量以后，便可以在运算中使用 In4：f2 Out4=32In5：f1+bOut5=(1+b)2In6：gz,4 Out6=5 如果忘记了已定义的函数的内容，可以使用？f查询f的定义当函数或变量使用完以后，最好将其清除，以免带来麻烦 3符号运算 符号运算即代数式的运算它是Mathemaica的重要功能下面简介符号运算的主要功能 （1）符号赋值 Mathemaica不仅可以把一个常值赋给一个符号，还可以把一个表达式赋给一个符号其规则如下： xvalue-将value赋给x x-清除赋给x的值 expr/.x- value -用value 替换expr中的x expr/.x-xvalue,y-yvalue-用xvalue,yvalue分别替换expr中的x,y. 例： In1：tlx Out11+x In2： l- t 2 Out21-(1+x)2 In3：t Out31-(1+x)2 In4：l- t 2 Out41-t2 In5：%2/.x-2Out5-8（2）代数式变换 Mathernatica提供了许多进行代数式变换的一些函数，下面列出常用的函数 Expandexpr-展开expr ExpandAllexpr-展开expr的分子、分母 Factorexpr-对expr进行因式分解 Togetherexpr-对expr进行通分 Apartexpr-将 expr分解为简单分式 Cancelexpr-消去exp r的分子、分母的公因式 Simplifyexpr-把expr化为最少项形式 例： In1：t=(x-1)2(2+x)/(1+x)(x-3)2)， In2：Expandt （展开分子，分母不变） In3：ExpandAllt （展开分子、分母） In4：=Together （通分） In5：Apart （化为部分分式） In6：Factor% （分解因式） In7：Simplify%5 (将表达式化简) 除了上述常用的变换外，Mathematica还可以进行许多种类型的变换下面再看一些例子In8：Expand2Cosx3*Sin2x2, Trig-True （展开三角函数） Out8：In9：Factor%,Trig-TrueOut9=8 Cosx5Sinx2In10:=ComplexExpandSinx+y*I （展开复函数） Out10：CoshySinx+ICosxSinyIn11：s=Expand(x+y)3；In12：Coefficients,x2 (取出s中x2项的系数)Out12=3y In13：=Numerator%1 (取出%1中的分子)Out13=(-1+x)2(2+x)In14：=Denominator%1 (取出%1中的分母)Out14=(-3+x)2(1+x)Mathematica还允许用户自己定义变换规则，例如： In15：mysin=Sin2*x_-2SinxCosx； In16：Sin2*(x+y)2/.mysin Out16=2Cos(x+ y)2Sin(x+ y)2 总之Mathematica进行变换的功能是非常强的（3）解方程 Mathematica可以用多种方法求解符号方程下面列出主要的解法： Solveequ,vars-求方程的一般解 Reduceequ,vars-求方程的全部解 NSolveequ,vars-求方程的数值解 FindRootequ,x,a-求方程在 a附近的数值解其中，equ是待求解的方程，var是未知量例 In1:=Solvea*x+b=0,x注：方程中，等号必须用“= =”Out1=x-b/aIn2:=Reducea*x+b=0,xOut2=a = 0 & b = 0 | x =-b/a & a != 0使用Reduce给出了a!=0时的解和a=0,b=0时的解，（此时x为任意值） 对四次及四次以下的代数方程， Mathematica总能给精确解四次以上的方程，若能分解因式，亦可给出精确解 In3：=Solvex3+3x2+ 3x+ 2= 0，x Out3= 当求不出精确解时，Mathemaica以符号形式给出结果 In4：=x5+5x+10； In5：=Solve4，x Out5=上述方程求不出精确解，此时可求数值解In6：=NSolve4，xOut6= 如果要求在某点附近的数值解，使用FindRoot In7：FindRootx*Sinx=1/2,x,1 Out7x-0.740841使用 Solve还可以求解方程组 1n8：=Solvex2+y2=1,x+y=a,x,yOut8三 微积分进行高等数学中的各种运算是Mathematica的主要功能Mathematica可进行微积分、线性代数和工程数学中的许多运算特别是其符号运算能力，令人惊叹现在Mathematica已受到越来越多科技工作者的欢迎和使用。1极限、微分和积分 微积分等主要运算：例 In1:=DSinx2,xOut1=2xCosx2In2:=Dxn,x,3Out2= In3:=Dy3*Logx+y,x,yOut3= 也可以求出抽象函数的导数In4:=Dx*fx5,xOut4= 求不定积分，对Mathematica而言易如反掌In5:=Integrate1/(x4-1),xOut5= 可以验证In6:=SimplifyD%,xOut6= 求定积分In7:= IntegrateLogx,x,a,bOut7=a-b-aLoga+bLogb 也可以通过File-Palettes-BasicCalculations输入In8:= Out8=a-b-aLoga+bLogbIn9:= Integratex*x+y*y,x,0,1,y,0,Sqrt1-x*xOut9=2. 函数的幂级数展开Mathematica作幂级数展开可达到任意精度。进行幂级数展开，使用下述函数：Seriesexpr,x,x0,n- expr在x=x0点的n阶幂级数展开式Normalseries-去掉展开式的余项例In1:=SeriesLog1+x,x,0,5Out1= In2:= Normal%Out2=抽象函数展开In3:=Seriesfx,x,0,4Out3=3微分方程求解微分方程与代数方程类似,导数用单引号表示DSolveequ,yx,x-解yx的微分方程DSolveequ,xt,yt,t-解自变量为t的微分方程组例 In1:= DSolveyx=a*yx, yx,xOut1=如给出初始条件，就明确了待定常数In2:= DSolveyx=a*yx,y0=2, yx,xOut2= In3:= DSolvey”x+yx=1,y0=2,y0=5, yx,xOut3= In4:= DSolveyt=xt,xt=yt,x0=1,y0=2, xt,yt,tOut4= 四 线性代数1构造矩阵a11,a1n,a21,a2n,am1,amn-构造mn矩阵Tablefi,j,i,m,j,n-构造mn矩阵Arrayf,m,n-以fi,j为元素，构造mn矩阵DiagonalMatrixa1,a2,an-构造n阶对角矩阵IdentityMatrixn -构造n阶单位矩阵TableIfi=j,fi,j,0,i,m,j,n-构造下三角矩阵mi,j-取矩阵m的元素m(i,j)mi-取矩阵m的第i行Transposemk -取矩阵m的第k列例 In1:=a,b,c,d Out1= a,b,c,d如果希望得到矩阵形式，可使用函数MatrixFormIn2:=MatrixForm% Out2=/ MatrixForm=In3:=Table3*(i-1)+j,i,3,j,3Out3=1,2,3,4,5,6,7,8,9In4:= MatrixFormDiagonalMatrixa,b,cOut4=/ MatrixForm=In5:=p=Arraya,2,3Out5= 取出第2行In6:=p2Out6=a2,1,a2,2,a2,3取出第3列In7:= Transposep3Out7=a1,3,a2,32矩阵运算Mathematica可进行矩阵的各种运算,如矩阵求逆、矩阵的转置、矩阵与向量的乘法等.下面列出主要的运算.记k为常数,u,v为向量,A,B为矩阵k*A-常数乘矩阵k+u-向量u的每一个元素加上ku+v-向量的对应元素相加u.v-向量的内积u*v-向量的对应元素相乘A.u-矩阵乘向量u.A-向量乘矩阵A.B-矩阵乘矩阵MatrixPowerA,n -矩阵乘幂TransposeA-求矩阵A的转置阵InverseA-求矩阵A的逆矩阵DetA-求矩阵A的行列式EigenvaluesA-求数字阵A的特征值EigentvectorsA-求数字阵A的特征向量LinearSolveA,v-求解线性方程组Ax=vChop%n-舍去第n个输出中无实际意义小量矩阵可以左乘以向量或右乘以向量, Mathematica也不区分“行”,或“列”向量,自动进行可能的运算.例:In1:=A=a,b,c,d; v=x,y;In2:=A.v (A左乘以v)Out2=ax+by,cx+dyIn3:=v.A (A右乘以v)Out3=ax+cy,bx+dyIn4:=InverseAOut4= 如果矩阵的元素是近似数，则求出的逆矩阵也是近似的。In5:=B=1.2,5.7,4.2,5.6; InverseBOut5=In6:=%.BOut6= 结果与单位矩阵有微小误差，用函数Chop消去无实际意义小量In7:=Chop%Out7=1.,0,0,1.前面已介绍了用Solve解线性方程组，但对于矩阵形式Ax=v的线性方程组，用LinearSolveA,v更方便.In8:=M=2,1,1,4; LinearSolveM,a,bOut8=五 数据拟合与优化1数据拟合Fitdata,funs,vars-用变量为vars,函数为funs，按最小二乘法拟合一组数据dataInterpolatingPolynomialx1,f1,x2,f2,x-多项式插值例In1:=d1=TableExp-x/2.,x,8 (生成一个数据表)Out1= In2:=Fitd1,1,x,x2,x3,x (用三次多项式拟合)Out2=In3:= Fitd1,Sinx,Exp-x/2,x （用三角函数与指数函数拟合）Out3= In4:=Chop%Out4= In5:=d2=FlattenTablex,y,1-2x+x*y+y2,x,2,5,y,10,15,1 (生成一个二元数据表，即一串三维点)Out5= In6:=Fitd2,1,x,y,x*y,y2,x,y （二元函数拟合）Out6= In7:=Chop%Out7=In8:= d5=0,2.5,0.5,2.1,1.,2.8,1.5,3.2;InterpolatingPolynomiald5,xOut8= 2优化FindMinimumf,x,x0-从x=x0始，求一元函数f的局部极小值FindMinimumf,x,x0,x1-从x0,x1始，求一元函数f的局部极小值FindMinimumf,x,x0,y,y0,- 从点（x0，y0,）始，求多元函数f的局部极小值ConstrainedMinf,inequalities,x,y,，求线性不等式约束下的f最小值,假定变量非负ConstrainedMaxf,inequalities,x,y,，求线性不等式约束下的f最大值,假定变量非负LinearProgrammingc,A,b-求解线性规划mincx|Ax=b,x=0例In1:= FindMinimumx4+3x2*y+5y2+x+y,x,0.1,y,0.2Out1= In2:= ConstrainedMax17x-20y+18z,x-y+z=15,x5,x+yPointSize0.02,RGBColor0,0,1,DisplayFunctionIdentity;pt1=ListPlotdata1,PlotJoined-True,PlotStyle-RGBColor0,0,1,DisplayFunctionIdentity;data2=1.2,5.8,1.8,5.6,2.7,7.9,3.8,8,4.5,4;pd2=ListPlotdata2,PlotStyle-PointSize0.02,RGBColor0,1,0.5,DisplayFunctionIdentity;pt2=ListPlotdata2,PlotJoined-True,PlotStyle-RGBColor0,1,0.5,DisplayFunctionIdentity;Showpd1,pt1,pd2,pt2,DisplayFunction$DisplayFunctionIn2:= p1=PlotSin2x,Cos3x,x,0,Pi,PlotStyleRGBColor1,0,0In3:=p2=ParametricPlotCost3,Sint3,t,0,2PiIn4:=Showp1,p2In5:=ListContourPlot22,33,55,32,42,23,34,44,52,43,52,452立体图形ListPlot3DArray- 画三元数组散点图Plot3Df,x, xmin,xmax , y, ymin,ymax -画二元函数曲面图ParametricPlot3Dfx,fy,fz,t,min,max-画参数方程表示的空间曲线ParametricPlot3Dfx,fy,fz,t,t1,t2,u,u1,u2-画参数方程表示的曲面图例 In1:= da=3,4,2,5,6,8,7,9,4;ListPlot3DdaIn2:=Plot3DExp-(x2+y2),x,-2,2,y,-2,2,Axes-False,Boxed-FalseIn3:=Plot3Dx2-y2,x,-3,3,y,-3,3,Axes-False In4:=ParametricPlot3DSint,Cost,t/8,t,0,15,Axes-FalseIn5:=ParametricPlot3DCost(3+Cosu),Sint(3+Cosu),Sinu,t,0,2Pi,u,0,2PiIn6:=ParametricPlot3D10SintCosu,10SintSinu,10Cost,t,0,Pi,u,0,2Pi 3声音Playf,t,min,max-根据函数f的值来播放声音In1:=PlaySin700t+250t*Sin15t,t,0,3七 编程1条件控制(1) Ifcond,then,else-如果cond成立，则执行then,否则执行else(2) Whichcond1,valu1, cond2,valu2,- 如果condi成立，则执行valuei(3) Switchexpr,form1,value1, form2,value2,-比较expr与每个formi,给出第一个匹配的value例 定义一个分段函数In1:=fx_:=Ifx=0 & x=1 & x0, n=Floorn/2;Printn 7313简单实例(1) 旋转曲面fx_ = Cosx + 1;Fori = 1, i -2, 2, -2, 2, -3, 3(2) 柱面y = Cosx; g1 = ParametricPlot3Dx, y, z, z, -0.01, 0, x, -Pi, PiFori = 1, i -Pi, Pi, -1, 1, 0, 5, DisplayFunction - Identity; Showg1, g2, DisplayFunction - $DisplayFunction, Boxed - False, Axes - None(3) 螺旋线rr = MovieParametricPlots Cos2 Pi s + t, s Sin2 Pi s + t, s, 0, 4, t, 0, 2Pi, Frames - 10, Axes - False, AspectRatio - Automatic, PlotRange - -4, 4, -4, 4; SpinShowrr, Frames - 10, SpinRange - 0 Degree, 180 Degree(4) 二分法fx_ = Sinx - x + 1; aa = 0.; bb = 10.; eps = 0.01;g1 = Plotfx, x, aa, bb, DisplayFunction - Identity;Fori = 1, i $DisplayFunction, PlotRange - 0, 10, f0, f10; Iffcc*faa 0, bb = cc, aa = cc(5) 正态分布mu = 0;normalfx_ = Exp-0.5(x - mu)/sigma)2/(Sqrt2Pisigma);Fork = 0, k -5, 5, 0, 0.4八 综合实例1最速降线O(0,0)A(x,y)B(1,1)xyC图2 质点沿曲线C下降一个在点O(0,0)静止的质点在重力作用下，该沿什么轨迹曲线C无摩擦地从点O(0,0)滑到点B(1,1),才能使所花时间T最短？(见图2)由于质点在下降时所增加的动能应等于所减少的势能，故质点在点A(x,y)处的速度v满足.其中，m为质点的质量。从而，且利用微积分的微元法可知T可由曲线积分求得： 从理论上推断，曲线C为摆线时T的值最小1，假如我们对此结论不太确信，则不妨对曲线C分别为摆线、直线、圆弧和抛物线的各种不同情况计算一下T的值并加以比较。（1）C为摆线。先让我们来考察当C为摆线时T的值。过点O与点B的摆线参数方程为其中，t=0与t=t0=2.412088分别对应点O与点B，用Mathematica软件求解，输入输出T=0.583203.（2）C为直线。我们总觉得，空间两点用直线连接时其路程最短，从而质点运动时间好象也应该是最短的。让我们来算一下此时的T值。过点O与点B的直线方程为x=y (0y1)，可见，的确是质点从点O到点B沿直线运动不如沿摆线运动快，T值增加了9.55%，还不算小呢。（3）C为圆弧。再让我们来考察当C为圆弧时T的值。此时过点O与点B的圆弧方程为其中，(a,b)为圆心坐标，r为半径且r1。因为圆心必在直线段OB的垂直平分线上，故推得a+b=1,进而得，此时T是r的函数，显然r越大，C越接近直线。具体的曲线形状如何，可以通过作图观察。在Mathematica软件中输入：b=(1-Sqrt-1+2r*r)/2;PlotT,r,1,3图3 T与r的关系可得图3,从图可见，T的最小值点在r=1.3附近。再输入：FindMinimumT,r,1.3，得输出：0.58512,r-1.33136这说明当曲线C为半径r=1.33136的圆弧时，质点从点O到点B的下降时间达到最小值0.58512，比C为直线时要快，但不如C为摆线时快。 （4）C为抛物线。再让我们来考察当C为抛物线时T的值。此时过点O与点B的抛物线方程为x=ay2+(1-a)y ，其中，参数a(0a1)影响抛物线的弯曲程度。，此时T是a的函数，显然a=0时，C为直线。我们可以作图观察T与a的关系。在Mathematica软件中输入：; PlotT,a,0,1图4 T与a的关系可得图4，可见T的最小值点在a=0.9附近。再输入：FindMinimumT,a,0.9，得输出：0.583778,a-0.913034这说明当曲线C为a=0.913034的抛物线时，质点从点O到点B的下降时间T达最小值0. 583778，比C为圆弧时要快，但不如C为摆线时快。但也已经十分接近，只是增加了0.1%. 以上从多个侧面验证了几种曲线都不如摆线的T值小，这使我们更加确信理论推导的结论。ABDC2追逐线一个追逐问题:如右图所示,正方形ABCD的四个顶点各有一人,分别是甲乙丙丁。在某一时刻(设为t0 =0),四人同时出发以匀速v走向顺时针方向的下一个人。如果他们始终保持对准目标,则最终将按螺旋状曲线汇合于中心点O.请求出这种情况下每个人的行进轨迹. 解 建立平面直角坐标系,取时间间隔为t,采样并计算每个人在每一时刻t的下一时刻t+t的位置(坐标)。设甲追逐乙, t时刻甲的坐标为(x1,y1), 乙的坐标为(x2,y2)。容易计算