无穷级数与微分方程(Mathematica)
119
项目四 无穷级数与微分方程
（仅供学习交流）
实验1 无穷级数（基础实验）
实验目的
观察无穷级数部分和的变化趋势,进一步理解级数的审敛法以及幂级数部分和对函数的 逼近. 掌握用Mathematica 求无穷级数的和, 求幂级数的收敛域, 展开函数为幂级数以及展 开周期函数为傅里叶级数的方法.
基本命令
1. 求无穷和的命令Sum
该命令可用来求无穷和. 例如,输入 Sum[1/n^2,{n,l,Infinity}]
则输出无穷级数的和为.6/2π 命令Sum 与数学中的求和号∑相当. 2. 将函数展开为幂级数的命令Series 该命令的基本格式为
Series[f[x],{x,x0,n}]
它将)(x f 展开成关于0x x -的幂级数. 幂级数的最高次幂为,)(0n x x -余项用10)(+-n x x 表 示. 例如,输入
Series[y[x],{x,0,5}] 则输出带皮亚诺余项的麦克劳林级数
[][][]()[]()[]()[][]6
5443320120
1024106102100x o x y x y x y x y x y y +++
+''+'+ 3. 去掉余项的命令Normal
在将)(x f 展开成幂级数后, 有时为了近似计算或作图, 需要把余项去掉. 只要使用 Normal 命令. 例如,输入
Series[Exp[x],{x,0,6}] Normal[%] 则输出
76
5432]x [O !
6x !5x !4x !3x !2x x 1++++++
+ !
6x 5x 4x !3x !2x x 16
5432++++++ 4. 强制求值的命令Evaluate
如果函数是用Normal 命令定义的, 则当对它进行作图或数值计算时, 可能会出现问题. 例如,输入
fx=Normal[Series[Exp[x],{x,0,3}]] Plot[fx,{x,-3,3}]
则只能输出去掉余项后的展开式 6
x 2x x 13
2+++
120 而得不到函数的图形. 这时要使用强制求值命令Evaluate, 改成输入 Plot[Evaluate[fx],{x,-3,3}] 则输出上述函数的图形.
5. 作散点图的命令ListPlot
ListPlot [ ]为平面内作散点图的命令, 其对象是数集,例如,输入
ListPlot[Table[j^2,{j,16}],PlotStyle->PointSize[0,012]]
则输出坐标为}16,16{,},3,3{},2,2{},1,1{2222 的散点图(图1.1).
图1.1
6. 符号“/；”用于定义某种规则,“/；”后面是条件. 例如,输入
Clear[g,gf];
g[x_]:=x/;0<=x<1 g[x_]:=-x/;-1<=x<0 g[x_]:=g[x –2]/;x>=1
则得到分段的周期函数
??
?
??≥-<≤<≤--=1x ),2x (g 1x 0,x 0x 1,
x )x (g
再输入
gf=Plot[g[x],{x,-1,6}] 则输出函数)(x g 的图形1.2.
图1.2
注:用Which命令也可以定义分段函数, 从这个例子中看到用“…(表达式)/; …(条件)”来定义周期性分段函数更方便些. 用Plot命令可以作出分段函数的图形, 但用Mathematica命令求分段函数的导数或积分时往往会有问题. 用Which定义的分段函数可以求导但不能积分. Mathematica内部函数中有一些也是分段函数. 如:Mod[x,1],Abs[x],Floor[x]和UnitStep[x]. 其中只有单位阶跃函数UnitStep[x]可以用Mathematica命令来求导和求定积分. 因此在求分段函数的傅里叶系数时, 对分段函数的积分往往要分区来积. 在被积函数可以用单位阶跃函数UnitStep的四则运算和复合运算表达时, 计算傅里叶系数就比较方便了.
实验举例
数项级数
例1.1 (教材例1.1)
(1)观察级数∑∞
=1
2 1
n n
的部分和序列的变化趋势.
(2) 观察级数∑∞
=11
n n
的部分和序列的变化趋势.
输入
s[n_]=Sum[1/k^2,{k,n}];data=Table[s[n],{n,100}];
ListPlot[data];
N[Sum[1/k^2,{k,Infinity}]]
N[Sum[1/k^2,{k,Infinity}],40]
则输出(1)
级数的近似值为1.64493.
输入
s[n_]=Sum[1/k,{k,n}];data=Table[s[n],{n,50}];
ListPlot[data,PlotStyle->PointSize[0.02]];
则输出(2)中级数部分和的的变化趋势图1.4.
121
122
例1.2 (教材 例1.2) 画出级数
∑
∞
=--1
1
1
)1(n n n
的部分和分布图. 输入命令
Clear[sn,g];sn=0;n=1;g={};m=3;
While[1/n>10^-m,sn=sn+(-1)^(n-1)/n;
g=Append[g,Graphics[{RGBColor[Abs[Sin[n]],0,1/n],
Line[{{sn,0},{sn,1}}]}]];n++];
Show[g,PlotRange->{-0.2,1.3},Axes->True];
则输出所给级数部分和的图形（图1.5）,从图中可观察到它收敛于0.693附近的一个数.
例1.3 求
∑∞
=++1
2
3
841
n n n
的值.
输入
Sum[x^(3k),{k,1,Infinity}] 得到和函数
3
3
x 1x +--
例1.4 (教材 例1.3) 设,!
10n a n
n = 求∑∞
=1
n n
a
.
输入
Clear[a];
a[n_]=10^n/(n!);
123
vals=Table[a[n],{n,1,25}];
ListPlot[vals,PlotStyle->PointSize[0.012]]
则输出n a 的散点图（1.6）,从图中可观察n a 的变化趋势. 输入 Sum[a[n],{n,l,Infinity}]
图1.6
求幂级数的收敛域 例1.5 (教材 例1.4) 求∑
∞
=+-0
21
)3(4n n
n n x 的收敛域与和函数.
输入
Clear[a];
a[n_]=4^(2n)*(x-3)^n/(n+1); stepone=a[n+1]/a[n]//Simplify
则输出
n
2)
x 3)(n 1(16++-+
再输入
steptwo=Limit[stepone,n->Infinity] 则输出
)x 3(16+-
这里对a[n+1]和a[n]都没有加绝对值. 因此上式的绝对值小于1时, 幂级数收敛; 大于1 时发散. 为了求出收敛区间的端点, 输入
ydd=Solve[steptwo==1,x] zdd=Solve[steptwo==-1,x]
则输出
????????????→??????????
??→1647x 1649x 与
由此可知,当
16
491647<
为了判断端点的敛散性, 输入 Simplify[a[n]/.x->(49/16)] 则输出右端点处幂级数的一般项为
124 1
n 1
+ 因此,在端点16
49
=x 处,级数发散. 再输入
Simplify[a[n]/.x->(47/16)] 则输出左端点处幂级数的一般项为
1
n )1(n
+- 因此,在端点16
47
=x 处, 级数收敛.
也可以在收敛域内求得这个级数的和函数. 输入 Sum[4^(2n)*(x-3)^n/(n+1),{n,0,Infinity}] 则输出
)
x 3(16)]
x 3(161[Log +-+---
函数的幂级数展开
例1.6 (教材 例1.5) 求x cos 的6阶麦克劳林展开式. 输入
Series[Cos[x],{x,0,6}] 则输出
76
42]x [o 720
x 24x 2x 1+-+-
注:这是带皮亚诺余项的麦克劳林展开式.
例1.6 (教材 例1.6) 求x ln 在1=x 处的6阶泰勒展开式.
输入
Series[Log[x],{x,1,6}] 则输出
.]x [o 6
)1x (5)1x (4)1x (3)1x (2)1x ()1x (76
5432+---+---+---
例1.7 (教材 例1.7) 求x arctan 的5阶泰勒展开式. 输入
serl=Series[ArcTan[x],{x,0,5}]; Poly=Normal[serl]
则输出x arctan 的近似多项式
5
x 3x x 5
3+-
通过作图把x arctan 和它的近似多项式进行比较. 输入
Plot[Evaluate[{ArcTan[x],Poly}],{x,-3/2,3/2},
PlotStyle->{Dashing[{0.01}],GrayLevel[0]},AspectRatio->l]
则输出所作图形（图1.7）, 图中虚线为函数x arctan ,实线为它的近似多项式.
125
例1.9 求()()2
211+--x x e 在1=x 处的8阶泰勒展开, 并通过作图比较函数和它的近似多项式. 输入
Clear[f];
f[x_]=Exp[-(x-1)^2*(x+1)^2];
poly2=Normal[Series[f[x],{x,1,8}]]
Plot[Evaluate[{f[x],poly2}],{x,-1.75,1.75},PlotRange->
{-2,3/2},PlotStyle->{Dashing[{0.01}], GrayLevel[0]}]
则得到近似多项式和它们的图1.8.
()()()()++-++-++--+--5
4
3
2
1161714141x x x x
()()()8761
17312814
x x x +--+--+- 图1.8
126 例 1.10 求函数x sin 在0=x 处的91,,7,5,3 阶泰勒展开, 通过作图比较函数和它的近似多项式, 并形成动画进一步观察. 因为
()())(!121sin 221
20
++=++-=
∑
n k n
k k
x o k x x 所以输入
Do[Plot[{Sum[(-1)^j*x^(2j+1)/(2j+1)!,{j,0,k}],
Sin[x]},{x,-40,40},
PlotStyle->{RGBColor[1,0,0],RGBColor[0,0,1]}],{k,1,45}]
则输出为x sin 的3阶和91阶泰勒展开的图形. 选中其中一幅图形,双击后形成动画. 图1.9是最后一幅图.
例1.11 利用幂级数展开式计算5240(精确到1010-). 因为
.311332432405
1455
??
? ??
-=-=
根据m x )1(+在0=x 处的展开式有
??
?
??-???-???-?-= 12
38245
31!3594131!2541315113240 故前)2(>n n 项部分和为
??????
? ?
???--?-=∑
∏
-===12411
431!
515315113n k k k
k i n k i S 输入命令
s[n_]=3(1-1/(5*3^4)-Sum[Product[5i-1,{i,1,k-1}]/(5^k k!3^(4k)),{k,2,n-1}]); r[n_]=Product[5i-1,{i,1,n-1}]/5^n/n!3^(4n-5)/80;
127
delta=10^(-10);n0=100;
Do[Print["n=",n,",","s[n]=",N[s[n],20]];
If[r[n] 则输出结果为 .992555739.22405 ≈ 傅里叶级数 例1.12 (教材 例1.8) 设)(x g 是以π2为周期的周期函数,它在[]ππ,-的表达式是 ? ??<≤<≤--=ππx x x g 0,10 ,1)( 将)(x g 展开成傅里叶级数. 输入 Clear[g]; g[x_]:=-1/;-Pi<=x<0 g[x_]:=1/;0<=x Plot[g[x],{x,-Pi,5 Pi},PlotStyle->{RGBColor[0,1,0]}]; 则输出)(x g 的图形 (图1.10). 因为)(x g 是奇函数, 所以它的傅里叶展开式中只含正弦项. 输入 b2[n_]:=b2[n]=2 Integrate[1*Sin[n*x],{x,0,Pi}]/Pi; fourier2[n_,x_]:=Sum[b2[k]*Sin[k*x],{k,1,n}]; tu[n_]:=Plot[{g[x],Evaluate[fourier2[n,x]]}, {x,-Pi,5 Pi}, PlotStyle->{RGBColor[0,1,0],RGBColor[1,0.3,0.5]},DisplayFunction->Identity]; (*tu[n]是以n 为参数的作图命令*) tu2=Table[tu[n],{n,1,30,5}]; (*tu2是用Table 命令作出的6个图形的集合*) toshow=Partition[tu2,2]; (*Partition 是对集合tu2作分割, 2为分割的参数*) Show[GraphicsArray[toshow]] (*GraphicsArray 是把图形排列的命令*) 则输出6个排列着的图形（图1.11）,每两个图形排成一行. 可以看到n 越大, )(x g 的傅里叶级数的前 n 项和与)(x g 越接近. 128 实验习题 1.求下列级数的和: (1) ;21 ∑ ∞ =k k k (2);) 12(1 1 2 ∑∞ =-k k (3) ;) 2(1 1 2 ∑∞ =k k (4) .)1(1 1 ∑ ∞ =--k k k 2. 求幂级数 ∑ ∞ =+--0 1 2)5() 1(n n n x 的收敛域与和函数. 3. 求函数)1ln()1(x x ++的6阶麦克劳林多项式. 4. 求x arcsin 的6阶麦克劳林多项式. 5. 设1 )(2+=x x x f ,求)(x f 的5阶和10阶麦克劳林多项式,把两个近似多项式和函数的 图形作在一个坐标系内. 6. 设)(x f 在一个周期内的表达式为??? ??<≤--=2121 1)(2x x x f , 将它展开为傅里叶级 数(取6项), 并作图. 7. 设)(x f 在一个周期内的表达式为? ??<≤-<≤=21,210, 1)(x x x x f , 将它展开为傅里叶级数 (取8项), 并作图. 8. 求级数∑ ∞ =1 sin k k k 的和的近似值. 129 实验2 微分方程（基础实验） 实验目的 理解常微分方程解的概念以及积分曲线和方向场的概念,掌握利用 Mathematica 求微分方程及方程组解的常用命令和方法. 基本命令 1. 求微分方程的解的命令DSolve 对于可以用积分方法求解的微分方程和微分方程组,可用Dsolve 命令来求其通解或特解. 例如,求方程023=+'+''y y y 的通解, 输入 DSolve[y ''[x]+3y '[x]+2y[x]==0,y[x],x] 则输出含有两个任意常数C[1]和C[2]的通解: {}{} ]2[C e ]1[C e ]x [y x x 2--+→ 注:在上述命令中,一阶导数符号 ' 是通过键盘上的单引号 ' 输入的,二阶导数符号 '' 要 输入两个单引号,而不能输入一个双引号. 又如,求解微分方程的初值问题: ,10,6,03400='==+'+''==x x y y y y y 输入 Dsolve[{y''[x]+4 y'[x]+3y[x]==0,y[0]==6, y'[0]==10},y[x],x] (*大括号把方程和初始条件放在一起*) 则输出 {}{} x 2x 3e 148(e ]x [y +-→- 2. 求微分方程的数值解的命令NDSolve 对于不可以用积分方法求解的微分方程初值问题,可以用NDSolve 命令来求其特解.例如 要求方程 5.0,032=+='=x y x y y 的近似解)5.10(≤≤x , 输入 NDSolve[{y'[x]==y[x]^2+x^3,y[0]==0.5},y[x],{x,0,1.5}] (*命令中的{x,0,1.5}表示相应的区间*) 则输出 {{y->InterpolatingFunction[{{0.,1.5}},< >]}} 注:因为NDSolve 命令得到的输出是解)(x y y =的近似值. 首先在区间[0,1.5]内插入一系 列点n x x x ,,,21 , 计算出在这些点上函数的近似值n y y y ,,,21 , 再通过插值方法得到 )(x y y =在区间上的近似解. 3. 一阶微分方程的方向场 一般地,我们可把一阶微分方程写为 ),(y x f y =' 的形式,其中),(y x f 是已知函数. 上述微分方程表明:未知函数y 在点x 处的斜率等于函数 f 在点),(y x 处的函数值. 因此,可在Oxy 平面上的每一点, 作出过该点的以),(y x f 为斜率 的一条很短的直线(即是未知函数y 的切线). 这样得到的一个图形就是微分方程),(y x f y =' 的方向场. 为了便于观察, 实际上只要在Oxy 平面上取适当多的点,作出在这些点的函数的 切线. 顺着斜率的走向画出符合初始条件的解,就可以得到方程),(y x f y ='的近似的积分曲 线. 130 例如, 画出0)0(,12=-=y y dx dy 的方向场. 输入 < g1=PlotVectorField[{1,1-y^2},{x,-3,3},{y,-2,2}, Frame->True, ScaleFunction->(1&),ScaleFactor->0.16, HeadLength->0.01,PlotPoints->{20,25}]; 则输出方向场的图形（图2.1）, 从图中可以观察到, 当初始条件为2/10=y 时, 这个微分方程的解介于1-和1之间, 且当x 趋向于-∞或∞时, )(x y 分别趋向于1-与1. 下面求解这个微分方程, 并在同一坐标系中画出方程的解与方向场的图解. 输入 sol=DSolve[{y'[x]==1-y[x]^2,y[0]==0},y[x],x]; g2=Plot[sol[[1,1,2]],{x,-3,3},PlotStyle->{Hue[0.1],Thickness[0.005]}]; Show[g2,g1,Axes->None,Frame->True]; 则输出微分方程的解x x e e x y 2211)(++-=,以及解曲线与方向场的图形（图2.2）. 从图中可以看 到, 微分方程的解与方向场的箭头方向相吻合. 实验内容 用Dsolve 命令求解微分方程 例2.1 (教材 例2.1) 求微分方程 2 2x xe xy y -=+'的通解. 输入 Clear[x,y]; DSolve[y '[x]+2x*y[x]==x*Exp[-x^2],y[x],x] 或 DSolve[D[y[x],x]+2x*y[x]==x*Exp[-x^2],y[x],x] 则输出微分方程的通解: ??? ?????????+→--]1[C e x e 21]x [y 22x 2x 其中C[1]是任意常数. 例2.2 (教材 例2.2) 求微分方程0=-+'x e y y x 在初始条件e y x 21 ==下的特解. 输入 Clear[x,y]; DSolve[{x*y ' [x]+y[x]-Exp[x]==0,y[1]==2 E},y[x],x] 则输出所求特解: 131 ?? ? ??????????????? ??+→ x e e ]x [y x 例2.3 (教材 例2.3) 求微分方程x e y y y x 2cos 52=+'-''的通解. 输入 DSolve[y ''[x]-2y '[x]+5y[x]==Exp[x]*Cos[2 x],y[x],x]//Simplify 则输出所求通解: ??? ?????????-++→])x 2[Sin ])1[c 4x (2]x 2[Cos ])2[c 81((e 81]x [y x 例2.4 (教材 例2.4) 求解微分方程x e x y +=''2, 并作出其积分曲线. 输入 g1=Table[Plot[E^x+x^3/3+c1+x*c2,{x,-5,5}, DisplayFunction->Identity],{c1,-10,10,5},{c2,-5,5,5}]; Show[g1,DisplayFunction->$DisplayFunction]; 图2.3 例2.5 (教材 例2.5) 求微分方程组???????=--=++02y x dt dy e y x dt dx t 在初始条件0,100====t t y x 下的特 解. 输入 Clear[x,y,t]; DSolve[{x' [t]+x[t]+2 y[t]==Exp[t], y'[t] -x[t]- y[t]==0, x[0]==1,y[0]==0},{x[t],y[t]},t] 则输出所求特解: ??? ??????? ??+-→→])t [Sin ]t [Cos e (21]t [y ],t [Cos ]t [x t 例2.6 验证c y y x =+--)3305(1515 2是微分方程2 )(42-='y x x y 的通解. 输入命令 132 < sol=(-5x^3-30y+3y^5)/15==C; g1=ImplicitPlot[sol/.Table[{C->n},{n,-3,3}],{x,-3,3}]; g2=PlotVectorField[{1,x^2/(y^4-2)},{x,-3,3},{y,-3,3}, Frame->True,ScaleFunction->(1&),ScaleFactor->0.16, HeadLength->0.01,PlotPoints->{20,25}]; g=Show[g2,g1,Axes->None,Frame->True]; Show[GraphicsArray[{g1,g2,g}]]; 则分别输出积分曲线如图2.4(a), 微分方程的方向场如图2.4(b). 以及在同一坐标系中画出积分曲线和方向场的图形如下图2.4 (c). 图2.4 从图 2.4(c)中可以看出微分方程的积分曲线与方向场的箭头方向吻合, 且当∞→x 时, 无论初始条件是什么, 所有的解都趋向于一条直线方程. 例2.7 (教材 例2.6) 求解微分方程,)1(1 22/5+=+-x x y dx dy 并作出积分曲线. 输入 < DSolve[y' [x]-2y[x]/(x+1)==(x+1)^(5/2),y[x],x] 则输出所给积分方程的解为 ??? ??????? ??+++→]1[C )x 1()x 1(32]x [y 22/7 下面在同一坐标系中作出这个微分方程的方向场和积分曲线(设),3,2,1,0,1,2,3---=C 输入 t=Table[2(1+x)^(7/2)/3+(1+x)^2c,{c,-1,1}]; g1=Plot[Evaluate[t],{x,-1,1},PlotRange->{{-1,1},{-2,2}}, PlotStyle->RGBColor[1,0,0],DisplayFunction->Identity]; g2=PlotVectorField[{1,-2y/(x+1)+(x+1)^(5/2)},{x,-0.999,1},{y,-4,4}, Frame->True,ScaleFunction->(1&), ScaleFactor->0.16, HeadLength->0.01, PlotPoints->{20,25},DisplayFunction->Identity]; Show[g1,g2,Axes->None,Frame->True,DisplayFunction->$DisplayFunction]; 133 则输出积分曲线的图形（图2.5）. 图2.5 例2.8 求解微分方程,2)21(22-+='-y x y xy 并作出其积分曲线. 输入命令 < DSolve[1-2*x*y[x]*y' [x]==x^2+(y[x])^2-2,y[x],x] 则得到微分方程的解为 .)2(3 23C y x x y ++-+= 我们在33≤≤-C 时作出积分曲线, 输入命令 t1=Table[(3+Sqrt[3])Sqrt[3+24x^2-4x^4-4*c*x]/(6*x),{c,-3,3}]; t2=Table[(3-Sqrt[3])Sqrt[3+24x^2-4x^4-4*c*x]/(6*x),{c,-3,3}]; gg1=Plot[Evaluate[t1],{x,-3,3},PlotRange->{{-3,3},{-3,3}}, PlotStyle->RGBColor[1,0,0],DisplayFunction->Identity]; gg2=Plot[Evaluate[t2],{x,-3,3},PlotRange->{{-3,3},{-3,3}}, PlotStyle->RGBColor[1,0,0],DisplayFunction->Identity]; g1=ContourPlot[y-x^3/3-x*(-2+y^2),{x,-3,3},{y,-3,3},PlotRange->{-3,3},Contours->7,Cont ourShading->False,PlotPoints->50,DisplayFunction->Identity]; g2=PlotVectorField[{1,(x^2+y^2-2)/(1-2*x*y)},{x,-3,3},{y,-3,3}, Frame->True,ScaleFunction->(1&),ScaleFactor->0.16, HeadLength->0.01,PlotPoints->{20,25}, DisplayFunction->Identity]; Show[g1,g2,Axes->None,Frame->True, DisplayFunction->$DisplayFunction]; Show[gg1,gg2,g2,Axes->None,Frame->True, 134 DisplayFunction->$DisplayFunction]; 则输出微分方程的向量场与积分曲线, 并输出等值线的图2.6. 图2.6 用NDSolve 命令求微积分方程的近似解 例2.9 (教材 例2.7) 求初值问题:1,0)1()1(2 .1=='-++=x y y xy y xy 在区间[1.2,4]上的近 似解并作图. 输入 fl=NDSolve[{(1+x*y[x])*y[x]+(1-x*y[x])*y'[x]==0,y[1.2]==1},y,{x,1.2,4}] 则输出为数值近似解(插值函数)的形式: {{y->InterpolatingFunction[{{1.2,4.}},< >]}} 用Plot 命令可以把它的图形画出来.不过还需要先使用强制求值命令Evalu-ate, 输入 Plot[Evaluate[y[x]/.fl],{x,1.2,4}] 则输出近似解的图形（图2.7）. 图2.7 如果要求区间[1.2,4]内某一点的函数的近似值, 例如8 .1=x y ,只要输入 y[1.8]/.fl 则输出所求结果 {3.8341} 135 例2.10 (教材 例2.8) 求范德波尔(Van der Pel)方程 5.0,0,0)1(0 2-=' ==+'-+''==x x y y y y y y 在区间[0,20]上的近似解. 输入 Clear[x,y]; NDSolve[{y''[x]+(y[x]^2-1)*y'[x]+y[x]==0,y[0]==0,y'[0]==-0.5},y,{x,0,20}]; Plot[Evaluate[y[x]/.%],{x,0,20}] 可以观察到近似解的图形（图2.8）. 图2.8 ?? ?? ?==+-'1)1(0 1sin 2y x y x y x 的数值解, 并作出数值解的图形. 输入命令 < sol=NDSolve[{x*y'[x]-x^2*y[x]*Sin[x]+1==0,y[1]==1}, y[x],{x,1,4}]; f[x_]=Evaluate[y[x]/.sol]; g1=Plot[f[x],{x,1,4},PlotRange->All, DisplayFunction->Identity]; g2=PlotVectorField[{1,(x^2*y*Sin[x]-1)/x},{x,1,4},{y,-2,9}, Frame->True,ScaleFunction->(1&),ScaleFactor->0.16, HeadLength->0.01,PlotPoints->{20,25}, DisplayFunction->Identity]; 136 g=Show[g1,g2,Axes->None,Frame->True]; Show[GraphicsArray[{g1,g}],DisplayFunction-> $DisplayFunction]; 则输出所给微分方程的数值解及数值解的图2.9. 例2.11 (教材 例2.9) 求出初值问题 ?? ?? ?='==+'+''0)0(,1)0(cos sin 22y y x y x y y 的数值解, 并作出数值解的图形. 输入 NDSolve[{y''[x]+Sin[x]^2*y'[x]+y[x]==Cos[x]^2, y[0]==1,y'[0]==0},y[x],{x,0,10}] Plot[Evaluate[y[x]/.%],{x,0,10}]; 则输出所求微分方程的数值解及数值解的图形（图2.10）. 图2.10 例2.12 (教材 例2.10) 洛伦兹(Lorenz)方程组是由三个一阶微分方程组成的方程组.这三个方程看似简单, 也没有包含复杂的函数, 但它的解却很有趣和耐人寻味. 试求解洛伦兹方程组 137 ,0 )0(,4)0(,12)0()(4)()()() ()(45)()()()(16)(16)(?????? ?===-='-+-=' -='z y x t z t y t x t z t y t x t z t x t y t x t y t x 并画出解曲线的图形. 输入 Clear[eq,x,y,z] eq=Sequence[x'[t]==16*y[t]-16*x[t], y'[t]==-x[t]*z[t]-y[t]+45x[t],z'[t]==x[t]*y[t]-4z[t]]; sol1=NDSolve[{eq,x[0]==12,y[0]==4,z[0]==0}, {x[t],y[t],z[t]},{t,0,16},MaxSteps->10000]; g1=ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.sol1],{t,0,16}, PlotPoints->14400,Boxed->False,Axes->None]; 则输出所求数值解的图形（图2.11(a)）. 从图中可以看出洛伦兹微分方程组具有一个奇异吸引子, 这个吸引子紧紧地把解的图形“吸”在一起. 有趣的是, 无论把解的曲线画得多长, 这些曲线也不相交 . 图2.11 改变初值为,10)0(,10)0(,6)0(=-==z y x 输入 sol2=NDSolve[{eq,x[0]==6,y[0]==-10,z[0]==10}, {x[t],y[t],z[t]},{t,0,24},MaxSteps->10000]; g2=ParametricPlot3D[Evaluate[{x[t],y[t],z[t]}/.sol2],{t,0,24}, PlotPoints->14400,Boxed->False,Axes->None]; Show[GraphicsArray[{g1,g2}]]; 则输出所求数值解的图形（图2.11(b)）. 从图中可以看出奇异吸引子又出现了, 它把解“吸”在某个区域内, 使得所有的解好象是有规则地依某种模式缠绕. 138 实验习题 1. 求下列微分方程的通解: (1) ;0136=+'+''y y y (2) ();024=+''+y y y (3) ;2sin 52x e y y y x =+'-'' (4) .)1(963x e x y y y +=+'-'' 2. 求下列微分方程的特解: (1) ;15,0,029400='==+'+''==x x y y y y y (2) .1,1,02sin =' ==++''==π π x x y y x y y 3. 求微分方程0cos 2)1(2=-+'-x xy y x 在初始条件10 ==x y 下的特解.分别求精确解 和数值解)10(≤≤x 并作图. 4. 求微分方程组???? ???=--=++t t e y x dt dy e y x dt dx 235的通解. 5. 求微分方程组??????==+-==-+==4,081,0300t t y y x dt dy x y x dt dx 的特解. 6. 求欧拉方程组324x y y x y x =-'+''的通解. 7. 求方程5,0,011='==+'+''==x x y y y y x y 在区间[0,4]上的近似解. 实验3 抛射体的运动（续）（综合实验） 实验目的 通过微分方程建模和Mathematica 软件,在项目一实验5的基础上,进一步研 究在考虑空气阻力的情况下抛射体的运动. 问题 根据侦察,发现离我军大炮阵地水平距离10km 的前方有一敌军的坦克群正以每小 时50km 向我军阵地驶来,现欲发射炮弹摧毁敌军坦克群. 为在最短时间内有效摧毁敌军坦 克,要求每门大炮都能进行精射击,这样问题就可简化为单门大炮对移动坦克的精确射击 问题. 假设炮弹发射速度可控制在0.2km/s 至0.5km/s 之间,问应选择怎样的炮弹发射速度和 怎样的发射角度可以最有效摧毁敌军坦克. 说明 本节我们研究受到重力和空气阻力约束的抛射体的射程. 用))(),((t y t x 记抛射体 的位置, 其中x 轴是运动的水平方向, y 轴是垂直方向. 通过在0=y 的约束下最大化x , 可以 计算出使抛射体的射程最大的发射角. 假设0=t 时抛射体(炮弹)在原点(0,0)以与水平线夹角 为,α初始速度为0v 发射出去. 它受到的空气阻力为 .,?? ? ??-=-=dt dy dt dx k kv F r (3.1) 重力为 Mathematica函数大全--运算符及特殊符号一、运算符及特殊符号 Line1;执行Line，不显示结果 Line1,line2顺次执行Line1，2，并显示结果 ?name关于系统变量name的信息 ??name关于系统变量name的全部信息 !command执行Dos命令 n! N的阶乘 !!filename显示文件内容 a-b减 a*b或a b 乘 a/b除 a^b 乘方 base^^num以base为进位的数 lhs&&rhs且 lhs||rhs或 !lha非 ++,-- 自加1，自减1 +=,-=,*=,/= 同C语言 >,<,>=,<=,==,!=逻辑判断（同c） lhs=rhs立即赋值 lhs:=rhs建立动态赋值 lhs:>rhs建立替换规则 expr//funname相当于filename[expr] expr/.rule将规则rule应用于expr expr//.rule 将规则rule不断应用于expr知道不变为止param_ 名为param的一个任意表达式（形式变量）param__名为param的任意多个任意表达式（形式变量） 二、系统常数 Pi 3.1415....的无限精度数值 E 2.17828...的无限精度数值 Catalan 0.915966..卡塔兰常数 EulerGamma 0.5772....高斯常数 GoldenRatio 1.61803...黄金分割数 Degree Pi/180角度弧度换算 I复数单位 Infinity无穷大 二阶线性常微分方程的幂级数解法 从微分方程学中知道，在满足某些条件下，可以用幂级数来表示一个函数。因此，自然想到，能否用幂级数来表示微分方程的解呢？ 例1、求方程 ''0y xy -=的通解 解：设2012n n y a a x a x a x =+++++…… 为方程的解，这里(0,1,2,,,)i a i n =……是待定常系数，将它对x 微分两次，有 ''212312132(1)(1)n n n n y a a x n n a x n na x --+=?+?++-+++ 将y ，'y 的表达式代入方程，并比较的同次幂的系数，得到 x -∞<<∞2210a ?=，30320,a a ?-= 41430,a a ?-= 52540,a a ?-= 或一般的可推得 32356(31)3k a a k k = ?????-? ， 1 3134673(31) k a a k k += ??????+ ， 320k a += 其中1a ，2a 是任意的，因而代入设的解中可得： 36347 01[1][] 2323562356(31)33434673(31) n x x x x x y a a x n n n n =+++++++++?????????-????????+ 这个幂级数的收敛半径是无限大的，因而级数的和（其中包括两个任意常数0a 及1a ）便是所要求的通解。 例6 求方程'''240y xy y --=的满足初值条件(0)0y =及'(0)1y =的解。 解 设级 2012n n y a a x a x a x =+++++……为方程的解。首先，利用初值 条件，可以得到 00a =， 11a =， 因而 2323'2123''223123232(1)n n n n n n y x a x a x a x y a x a x na x y a a x n n a x --=+++++=+++++=+?++-+ 将y ，'y ，''y 的表达式带入原方程，合并x 的各同次幂的项，并令各项系数等于零，得到 21422 0,1,0,,,1 n n a a a a a n -==== - 因而 567891111 ,0,,0,,2!63!4! a a a a a = ===== 最后得 21111 (1)!! k a k k k += ?=- , 20k a =, 对一切正整数k 成立。 将i a (0,1,2,)i = 的值代回2012n n y a a x a x a x =+++++……就得到 521 3 2!! k x x y x x k +=+++++ 2 422 (1),2!! k x x x x x xe k =++++ += 这就是方程的满足所给初值条件的解。 是否所有方程都能按以上方式求出其幂级数解？或者说究竟方程应该满足什么条件才能保证它的解可用幂级数来表示呢？级数的 Mathematica的内部常数 Pi , 或π(从基本输入工具栏输入, 或“Esc”+“p”+“Esc”）圆周率π E (从基本输入工具栏输入, 或“Esc”+“ee”+“Esc”）自然对数的底数e I (从基本输入工具栏输入, 或“Esc”+“ii”+“Esc”）虚数单位i Infinity, 或∞(从基本输入工具栏输入, 或“Esc”+“inf”+“Esc”）无穷大∞ Degree 或°(从基本输入工具栏输入,或“Esc”+“deg”+“Esc”）度 Mathematica的常用内部数学函数 指数函数Exp[x]以e为底数 对数函数Log[x]自然对数，即以e为底数的对数 Log[a，x]以a为底数的x的对数 开方函数Sqrt[x]表示x的算术平方根 绝对值函数Abs[x]表示x的绝对值 三角函数 （自变量的单位为弧度）Sin[x]正弦函数 Cos[x]余弦函数 Tan[x]正切函数 Cot[x]余切函数 Sec[x]正割函数 Csc[x]余割函数 反三角函数ArcSin[x]反正弦函数 ArcCos[x]反余弦函数 ArcTan[x]反正切函数 ArcCot[x]反余切函数 ArcSec[x]反正割函数 ArcCsc[x]反余割函数 双曲函数Sinh[x]双曲正弦函数 Cosh[x]双曲余弦函数 Tanh[x]双曲正切函数 Coth[x]双曲余切函数 Sech[x]双曲正割函数 Csch[x]双曲余割函数 反双曲函数ArcSinh[x]反双曲正弦函数 ArcCosh[x]反双曲余弦函数 ArcTanh[x]反双曲正切函数 ArcCoth[x]反双曲余切函数 ArcSech[x]反双曲正割函数 ArcCsch[x]反双曲余割函数 求角度函数ArcTan[x，y]以坐标原点为顶点，x轴正半轴为始边，从原点到点（x，y）的射线为终边的角，其单位为弧度 数论函数GCD[a,b,c,．．．]最大公约数函数 LCM[a,b,c,．．．]最小公倍数函数 Mathematica函数大全一、运算符及特殊符号 Line1; 执行Line，不显示结果 Line1,line2 顺次执行Line1，2，并显示结果 ?name 关于系统变量name的信息 ??name 关于系统变量name的全部信息 !command 执行Dos命令 n! N的阶乘 !!filename 显示文件内容 < >,<,>=,<=,==,!= 逻辑判断（同c） lhs=rhs 立即赋值 lhs:=rhs 建立动态赋值 lhs:>rhs 建立替换规则 lhs->rhs 建立替换规则 expr//funname 相当于filename[expr] expr/.rule 将规则rule应用于expr expr//.rule 将规则rule不断应用于expr知道不变为止 param_ 名为param的一个任意表达式（形式变量） param__ 名为param的任意多个任意表达式（形式变量） 二、系统常数 Pi 3.1415....的无限精度数值 E 2.17828...的无限精度数值 Catalan 0.915966..卡塔兰常数 EulerGamma 0.5772....高斯常数 GoldenRatio 1.61803...黄金分割数 Degree Pi/180角度弧度换算 I 复数单位 Infinity 无穷大 -Infinity 负无穷大 ComplexInfinity 复无穷大 Indeterminate 不定式 三、代数计算 Expand[expr] 展开表达式 Factor[expr] 展开表达式 Simplify[expr] 化简表达式 FullSimplify[expr] 将特殊函数等也进行化简 PowerExpand[expr] 展开所有的幂次形式 ComplexExpand[expr,{x1,x2...}] 按复数实部虚部展开 FunctionExpand[expr] 化简expr中的特殊函数 Collect[expr, x] 合并同次项 Collect[expr, {x1,x2,...}] 合并x1,x2,...的同次项 Together[expr] 通分 Apart[expr] 部分分式展开 Apart[expr, var] 对var的部分分式展开 Cancel[expr] 约分 ExpandAll[expr] 展开表达式 ExpandAll[expr, patt] 展开表达式 FactorTerms[poly] 提出共有的数字因子 FactorTerms[poly, x] 提出与x无关的数字因子 FactorTerms[poly, {x1,x2...}] 提出与xi无关的数字因子 Coefficient[expr, form] 多项式expr中form的系数 高阶线性微分方程常用解法简介 关键词：高阶线性微分方程 求解方法 在微分方程的理论中，线性微分方程是非常值得重视的一部分内容，这不仅 因为线性微分方程的一般理论已被研究的十分清楚，而且线性微分方程是研究非线性微分方程的基础，它在物理、力学和工程技术、自然科学中也有着广泛应用。下面对高阶线性微分方程解法做一些简单介绍. 讨论如下n 阶线性微分方程：1111()()()()n n n n n n d x d x dx a t a t a t x f t dt dt dt ---++++= （1），其中()i a t （i=1,2,3,,n ）及f(t)都是区间a t b ≤≤上的连续函数，如果 ()0f t ≡，则方程（1）变为 1111()()()0n n n n n n d x d x dx a t a t a t x dt dt dt ---++++= （2），称为n 阶齐次线性微分方程，而称一般方程（1）为n 阶非齐次线性微分方程，简称非齐次线性微分方程，并且把方程（2）叫做对应于方程（1）的齐次线性微分方程. 1.欧拉待定指数函数法 此方法又叫特征根法，用于求常系数齐次线性微分方程的基本解组。形如 111121[]0,(3),n n n n n n n d x d x dx L x a a a x dt dt dt ---≡++++=其中a a a 为常数，称为n 阶常系数齐次线性微分方程。 111111111111[]()()()n t n t t t t n n n n n n n t t n n n n n n n d e d e de L e a a a e dt dt dt a a a e F e F a a a n λλλλλλλλλλλλλλλλ---------≡++++=++++≡≡++++其中=0(4)是的次多项式. ()F λ为特征方程，它的根为特征根. 1.1特征根是单根的情形 设12,,,n λλλ是特征方程111()0n n n n F a a a λλλλ--≡++++=的n 个彼此不相等的根，则应相应地方程（3）有如下n 个解：12,,,.n t t t e e e λλλ（5）我们指出这n 个解在区间a t b ≤≤上线性无关，从而组成方程的基本解组. 如果(1,2,,)i i n λ=均为实数，则（5）是方程（3）的n 个线性无关的实值 解，而方程（3）的通解可表示为1212,n t t t n x c e c e c e λλλ=+++其中12,,,n c c c 为任意常数. 如果特征方程有复根，则因方程的系数是实常数，复根将称对共轭的出现.设1i λαβ=+是一特征根，则2i λαβ=-也是特征根，因而于这对共轭复根 表达式： Plot[4 x - 9, {x, 0, 9}] f[x_] = x^3 Plot[f[x], {x, 0, 9}] a = Plot[4 x - 9, {x, 0, 9}] b = Plot[x^3, {x, 0, 3}] 两图画在一个坐标系 Show[a, b] a = Plot[4 x - 9, {x, 0, 9}] b = Plot[x^3, {x, 0, 3}] 两图画在一起（一排） c = GraphicsArray[{a, b}] Show[c] a = Plot[4 x - 9, {x, 0, 9}] b = Plot[x^3, {x, 0, 3}] c = GraphicsArray[{a}, {b}] 两图画在一起（两排） Show[c] 二维画图： Automatic 默认值 DisplayFunction -> Identity 不出现图 DisplayFunction -> $DisplayFunction 出现图 PlotRange -> All 画出所有点，指定区域点 PlotStyle -> {RGBColor[1, 0, 0]} 图像颜色 PlotStyle -> {Dashing[{0.01}]} 图像成虚线 PlotStyle -> {Thickness[0.01]} 图像粗细 AxesLabel -> {"x/t", "y/cm"} 坐标标签 PlotLabel -> {"s-t"} 图像标签 Frame -> True 图像边框 Axes -> {True, True} 坐标轴的显示 AxesOrigin -> {0, -5} 设置坐标原点 GridLines -> {{-π, -π/2, 0, π/2, π}, {-1，-0.5，0, 0.5, 1}} 给坐标轴分网格 TextStyle -> {FontSize -> 30} 坐标字体大小AspectRatio -> Automatic 坐标比例一致 Ticks -> {{0, 1, 2, 3}, {0,10,20}} 在坐标轴上显示特定点ParametricPlot[x(t),y(t)},{t,0,6,}] 画参数方程 1 / 3 4.4 高阶微分方程降阶法、二阶线性微分方程幂级数解法 (Power series solution to second order linear ODE ) [教学内容] 1. 介绍高阶方程降阶法. 2. 介绍单摆方程及其椭圆积分函数.3. 介绍刘维尔公式求解二阶线性方程. [教学重难点] 重点是知道振幅反应(Amplitude Response ）； 难点是知道常见函数的拉普拉斯变换和逆变换. [教学方法] 预习1、2；讲授1、2 [考核目标] 1. 知道共振现象. 2. 知道拉普拉斯变换的概念和性质. 3. 知道常见函数的拉普拉斯变换和逆变换. 1. 高阶方程降阶法 例68. 数学摆方程及其求解 解：（1）模型描述：一根长度为l 的线一端是质量为m 的质点，另一端系于固定点O ，质点在垂直于地面的平面上作圆周运动。取逆时针运动方向作为摆与铅垂线所成角?的正方向， 质点运动加速度为22dt d m l ?，所受的力为?sin mg -. 于是单摆方程为??sin 22l g dt d -=. 下面考察如下柯西问题：??sin 22l g dt d -=，0)0(',)0(0==???. （2）令dt d v ?=，下面导出? d dv ，由??d dt dt dv d dv ?=知，dt d d dv dt dv dt d ???? ==22. 于是原方程化为 ??sin l g v d dv -=，这是一个一阶可分离变量型方程。 解得 C l g v +=?cos 212，再由初始条件0)0(',)0(0==???得到 )cos (cos 20??-± =l g v ，其中±号由摆运动位置确定. （3）将v 返回原变量得到 )cos (cos 20???-±=l g dt d ，这也是一个一阶可分离变量型方程。先考察摆从最大正角0?到0?-之间运动情形： )cos (cos 20???--=l g dt d l g t dt l g d t 22cos cos 000 -=-=-??? ? ???，特别地令?---=000 0cos cos 2????? d g l T ， 二阶线性常微分方程的幂级数解法 从微分方程学中知道，在满足某些条件下，可以用幂级数来表示一个函数。因此，自然想到，能否用幂级数来表示微分方程的解呢？ 例1、求方程 ''0y xy -=的通解 解：设2012n n y a a x a x a x =+++++…… 为方程的解，这里(0,1,2,,,)i a i n =……是待定常系数，将它对x 微分两次，有 ''212312132(1)(1)n n n n y a a x n n a x n na x --+=?+?+ +-+++ 将y ，'y 的表达式代入方程，并比较的同次幂的系数，得到 x -∞<<∞2210 a ?=，30320,a a ?-= 41430,a a ?-= 52540,a a ?-= 或一般的可推得 32356(31)3k a a k k = ?????-?， 1 3134673(31) k a a k k += ??????+， 320k a += 其中1a ，2a 是任意的，因而代入设的解中可得： 36 347 01[1][] 232356 2356(31)33434673(31) n x x x x x y a a x n n n n =+++ ++++++ ?????????-????????+ 这个幂级数的收敛半径是无限大的，因而级数的和（其中包括两个 任意常数0a 及1a ）便是所要求的通解。 例6 求方程'''240y xy y --=的满足初值条件(0)0y =及'(0)1y =的解。 解 设级数2012n n y a a x a x a x =+++++……为方程的解。首先，利用初值条件，可以得到 00a =， 11a =， 因而 2323'2123''223123232(1)n n n n n n y x a x a x a x y a x a x na x y a a x n n a x --=+++++=+++++=+?+ +-+ 将y ，'y ，''y 的表达式带入原方程，合并x 的各同次幂的项，并令各项系数等于零，得到 21422 0,1,0, ,,1 n n a a a a a n -==== - 因而 5678911 11,0,,0,,2!63!4! a a a a a = ===== 最后得 21111 (1)!! k a k k k += ?=- , 20k a =, 对一切正整数k 成立。 将i a (0,1,2,)i =的值代回2012n n y a a x a x a x =+++++……就得到 5 213 2! !k x x y x x k +=+++ ++ 2 4 22 (1),2! ! k x x x x x xe k =+++ ++= 这就是方程的满足所给初值条件的解。 是否所有方程都能按以上方式求出其幂级数解？或者说究竟方 第8章Mathematica中的常用函数8.1 运算符及特殊符号 Linel 执行Line，不显示结果 Linel,line2 顺次执行Line1，Line2，并显示结果 ?name 关于系统变量name的信息 ??name 关于系统变量name的全部信息 !command 执行Dos命令 N! N的阶乘 !!filename 显示文件内容 < -Infinity 负无穷大 Complexlnfinity 复无穷大 Indeterminate 不定式 8.3 代数计算 Expand[expr] 展开表达式 Factor[expr] 展开表达式 Simplify[expr] 化简表达式 FullSimplify[expr] 将特殊函数也进行化简PowerExpand[expr] 展开所有的幂次形式ComplexExpand[expr,{x1，x2…}] 按复数实部虚部展开FunctionExpand[expr] 化简表达式中的特殊函数 Collect[expr,x] 合并同次项 Collect[expr,{x1，x2，…}] 合并x1，x2,...的同次项 Together[expr] 通分 Apart[expr] 部分分式展开 Apart[expr,var] 对var的部分分式展开 Cancel[expr] 约分 ExpandAll[expr] 展开表达式 ExpandAll[expr,patt] 展开表达式 FactorTermsrpoly] 提出共有的数字因子 FactorTerms[poly，x] 提出与x无关的数字因子 FactorTerms[poly,(x1,x2…)] 提出与xi无关的数字因子 Coefficient[expr,form] 多项式expr中form的系数 Coefficient[expr,form，n] 多项式expr中form^n的系数 Exponent[expr,form] 表达式expr中form的最高指数 Numerator[expr] 表达式expr的分子 Denominator[expr] 表达式expr的分母 ExpandNumerator[expr] 展开expr的分子部分 8.4 解方程 Solve[eqns,vats] 从方程组eqns中解出Vats Solve[eqns,vats,elims] 从方程组eqns中削去变量elims，解出vats DSolve[eqn,y,x] 解微分方程，其中、y是x的函数 DSolve[{eqnl,eqn2,…}，{y1，y2…},] 解微分方程组，其中yi是x的函数DSolve[eqn,y,{x1,x2…}]解偏微分方程 Eliminate[eqns，Vats] 把方程组eqns中变量vars约去SolveAlways[eqns，vars] 给出等式成立的所有参数满足的条件Reduce[eqns，Vats] 化简并给出所有可能解的条件LogicalExpand[expr] 用&&和，，将逻辑表达式展开InverseFunction[f] 求函数f的反函数 Root[f,k] 求多项式函数的第k个根 Mathematica for Windows 常用用法 一、Mathematica 的主要功能 Mathematica 是美国Wolfram 公司开发的一个功能强大的计算机数学系统，提供了范围广泛的数学计算功能，主要包括三个方面：符号演算、数值计算、图形。例如：多项式的四则运算、展开、因式分解，有理式的各种计算，有理方程、超越方程的解，向量和矩阵的各种计算，求极限、导数、极值、不定积分、定积分、幂级数展开式，求解微分方程，作一元、二元函数的图形等等。 二、Mathematica 的基本知识 1．输入表达式：直接输入一个表达式（包括算式和命令，长表达式用“Enter ”换行）后，按“Shift+Enter ”执行，执行后以“Out[命令序号]= ……”形式输出执行结果，输出的结果可在后续的表达式中使用。 若命令后有分号，则不输出执行结果（图形输出与Print 命令除外）。 “%”表示上一个输出，“%%”表示倒数第2个输出，“%ｉ”表示第ｉ个 命令的输出。 2．运算符：+、-、*、/、^ ，“*”可用空格代替，“^”表示乘方。 如：In[1]:=2^10，输出为“Out[1]= 1024”，其中“In[1]:=”不需要输入。 In[2]:=3+5，Out[2]= 8；In[3]:=%-2，Out[3]= 6； In[4]:=%2+4，Out[4]= 12； In[5]:=1/3-1/4，Out[5]=12 1 ；In[6]:=N[%]，Out[6]= 0.0833333； In[7]:=N[%5+12,10]，Out[7]= 12.08333333（注意字母的大小写） 3．变量赋值：变量=表达式，“x=.”或Clear[x] 表示清除对x 的赋值。 表达式/.t ->c ，将表达式中的t 全替换成c 。?x ，查x 信息。 4．常用的数学常数：Pi (π)、E(e)、Infinity (∞)、I (1-) 5．常用的数学函数：Abs, Sin, Cos, Tan, Cot, ArcSin, Log (自然对数), Sqrt, Exp 如：In[1]:=Sqrt[2]+1；In[2]:=Sin[2]+ArcSin[1]；In[3]:=Exp[2]+% （自变量用［ ］括，区分大小写，首字母大写） 三、常用运算 1．多项式运算：In[1]:= (2+4*x^2)*(1-x)^3 或 In[1]:= t = (2+4*x^2)*(1-x)^3 （将右端表达式赋值给t ）； In[2]:=a=t/.x->4 (计算表达式t 当x=4时的值，并赋值给变量a ) In[3]:=a=. (清除变量a ) In[3]:=Expand[t]（展开）；In[4]:=Factor[%]（把上一个结果因式分解） 2．解方程：In[1]:=Solve[x^2+3*x = = 2]；In[2]:=N[%]； In[3]:=Solve[a*x-b= = 0, x]； In[4]:=NSolve[{x-2*y= =0,x^2-y= =1},{x,y}]（解方程组并得到数值解） 3．自定义函数：In[1]:= f [x_ ]:=x^2+2*x ； In[2]:=f[5]+7； In[3]:=f[a+b] 4．求极限：In[1]:=Limit[Sin[x]/x, x ->0]； In[2]:=Limit[(1+1/n)^n, n->Infinity]，Out[2]=E 5．求（偏）导数：In[1]:=D[a*x^2+3, x]；In[2]:=D[x^2+y^3-Sin[2*y], y]（对y 的偏导数）； In[3]:=D[Log[x], {x,2}] (求对x 的二阶导数)； In[4]:=D[Sin[x+y]*Exp[z*y^2],x,y] （求对x 、y 的二阶混合偏导数）； In[5]:=Simplify[%] (对前一结果化简)； In[6]:=D[Sin[x+y]*Exp[z*y^2],{x,2},{y,3}] 6．求不定积分：In[1]:=Integrate[x^2,x]；In[2]:=Integrate[1/(x^2+a^2),x] 7．定积分：In[1]:=Integrate[x^2, {x,0,1}]；In[2]:=Integrate[x^2,{x,a,b}]； In[3]:=Integrate[x^2+y^2, {x,0,a},{y,0,b}]；（求矩形域上的二重积分） In[4]:=Integrate[1, {x,-1,1},{y,-Sqrt[1-x^2],Sqrt[1-x^2]}]；Out[4]=Pi （圆面积） 8．幂级数展开：In[1]:=Series[Exp[x],{x,0,4}]（在x=0处展开到x 的四次幂） 9．矩阵的输入和输出：In[1]:= a ={{1,2},{3,4}}（定义一个2x2的矩阵a ，按 行写）； In[2]:=MatrixForm[a]（输出为矩阵形式）；In[3]:=Transpose[a]（a 的转置）； In[4]:=a[[2]]（a 的第2行）；In[5]:=Tanspose[a][[2]]（a 的第2列）； In[6]:=Inverse[a]（求a 的逆矩阵）；In[7]:=Det[a]（矩阵的行列式）； In[8]:=Eigenvalues[a]（求特征值）；In[9]:=Eigenvectors[a]（求特征向量）； In[10]:=RowReduce[a]（把a 化为阶梯形，可用于求矩阵的秩、判断线性相关性）； In[11]:= b ={{5,6,7},{8,9,10}}；In[12]:= a .b （矩阵a 与b 的乘积） 10．解线性方程组： In[1]:= a ={{3,4,5,6},{6,8,10,12},{4,5,6,7},{5,6,7,8}}；（a 的秩为2） In[2]:= b ={1,2,3,5}（列向量）；（增广矩阵的秩也为2） In[3]:=LinearSolve[a,b]（求线性方程组ax=b 的一个特解）； In[4]:=NullSpace[a]（求线性方程组ax=0的一个基础解系）； In[5]:= x =k1%4[[1]]+k2%4[[2]]+%3（ax=b 的全部解，k1、k2为任意常数） 11．求和：In[1]:=NSum[Sin[n]/n^3,{n,1,Infinity}]（求级数∑ ∞=13sin n n n 的和） 12．求极小值：In[1]:=FindMinimum[Sin[x]*Cos[x],{x,0.5}]（求函数在0.5附 近的极小值）； In[2]:=FindMinimum[Sin[x*y]*Exp[x^2],{x,0.2}, {y,0.3}]（求多元函数极小值） 13．求解线性规划问题：Min cx ，mx ≥b ，x ≥0，求向量x 。 In[1]:= c ={2,-3}（列向量）；In[2]:= m ={{-1,-1},{1,-1},{1,0}}； In[3]:= b ={-10,2,1}； In[4]:=LinearProgramming[c,m,b] 14．数据拟合：In[1]:= d ={{1,2.18},{1.2,2.56},{1.6,3.0},{1.8,2.66}}； In[2]:= f =Fit[d,{1, x, x^2}, x]（求和上面4个点吻合最好的二次多项式f ）； 检验效果：In[3]:=ListPlot[d]（画d 中4个点的图）； In[4]:=Plot[f,{x,0.8,2.0}]（画多项式f 在x 从0.8到2.0之间的图）； In[5]:=Show[%3, %4]（把上面两个图画在一起） 注：函数集{1, x, x^2}可以是更高次的或其它函数集，如三角函数集等。 15．一元函数作图：In[1]:=Plot[Exp[-x^2]*Sin[6*x],{x,-2,2}]（如图1） 参数方程作图：In[2]:=ParametricPlot[{Sin[t]^3,Cos[t]^3},{t,0,2*Pi}] 16．二元函数作图：In[1]:=Plot3D[Sin[x*y],{x,-Pi, Pi},{y,-Pi, Pi}]；（如图2） In[2]:=Plot3D[Sin[x*y],{x,-Pi, Pi},{y,-Pi, Pi},PlotPoints->40, ViewPoint->{2,-3,2}] In[3]:=ParametricPlot3D[{Cos[u]*Cos[v],Sin[u]*Cos[v],Sin[v]},{u,0,2*P i},{v,-Pi/2,Pi/2}] 17．数据画图：In[1]:= d ={{1,2},{3,4},{7,6}}；In[2]:=ListPlot[d]； In[3]:=ListPlot[d, PlotStyle->{RGBColor[1,0,0], PointSize[0.02]}]（红色 的大点）； 或直接用 In[4]:=ListPlot[{1,2},{3,4},{7,6}] 代替“In[2]:=”。 18．作图范围：In[1]:=Plot[x-x^3/6,{x,-4,4}]； In[2]:=Plot[x-x^3/6,{x,-4,4},PlotRange->{-5,2}]（限定纵坐标（函数值）范围） 19．图形组合：In[1]:=Plot[{Sin[x],Cos[x]},{x,0,2*Pi}]；或 In[2]:= g1=Plot[Sin[x],{x,0,2*Pi}, PlotStyle->{RGBColor[1,0,0]}]； In[3]:= g2=Plot[Cos[x],{x,0,2*Pi}, PlotStyle->{RGBColor[0,0,1]}]； In[4]:=Show[g1,g2]（把g1、g2画在一起） 20．文件的使用：In[1]:= y =25；In[2]:= a ={{1,4},{2,6}}；In[3]:= f [x_ ]:=x^2 ； In[4]:= g =Plot[Sin[x],{x,0,2*Pi}, PlotStyle->{RGBColor[1,0,0]}]； In[5]:=Save[“abc .m”,a,y,f,g]（将a, y, f, g 保存在文件“abc .m ”中，扩 展名为m ）； In[6]:=!!abc .m （显示文件内容）； In[1]:=< 高阶线性微分方程常用解法简介 摘要：本文主要介绍高阶线性微分方程求解方法，主要的内容有高阶线性微分方程求解的常 用方法如。 关键词：高阶线性微分方程 求解方法 在微分方程的理论中，线性微分方程是非常值得重视的一部分内容，这不仅 因为线性微分方程的一般理论已被研究的十分清楚，而且线性微分方程是研究非线性微分方程的基础，它在物理、力学和工程技术、自然科学中也有着广泛应用。下面对高阶线性微分方程解法做一些简单介绍. 讨论如下n 阶线性微分方程：1111()()()()n n n n n n d x d x dx a t a t a t x f t dt dt dt ---++++= （1），其中()i a t （i=1,2,3, ,n ）及f(t)都是区间a t b ≤≤上的连续函数，如果 ()0f t ≡，则方程（1）变为 1111()()()0n n n n n n d x d x dx a t a t a t x dt dt dt ---++++= （2），称为n 阶齐次线性微分方程，而称一般方程（1）为n 阶非齐次线性微分方程，简称非齐次线性微分方程，并且把方程（2）叫做对应于方程（1）的齐次线性微分方程. 1.欧拉待定指数函数法 此方法又叫特征根法，用于求常系数齐次线性微分方程的基本解组。形如 111121[]0,(3),n n n n n n n d x d x dx L x a a a x dt dt dt ---≡++++= 其中a a a 为常数，称为n 阶常系数齐次线性微分方程。 111111111111[]()()()n t n t t t t n n n n n n n t t n n n n n n n d e d e de L e a a a e dt dt dt a a a e F e F a a a n λλλλλλλλλλλλλλλλ---------≡++++=++++≡≡++++ 其中=0(4)是的次多项式. ()F λ为特征方程，它的根为特征根. 1.1特征根是单根的情形 设12,,,n λλλ 是特征方程111()0n n n n F a a a λλλλ--≡++++= 的n 个彼此不相等的根，则应相应地方程（3）有如下n 个解：12,,,.n t t t e e e λλλ （5）我们指出这n 个解在区间a t b ≤≤上线性无关，从而组成方程的基本解组. 如果(1,2,,)i i n λ= 均为实数，则（5）是方程（3）的n 个线性无关的实值解，而方程（3）的通解可表示为1212,n t t t n x c e c e c e λλλ=+++ 其中12,,,n c c c 为任意常数. 如果特征方程有复根，则因方程的系数是实常数，复根将称对共轭的出现. 微分方程的幂级数解法 函数是客观事物的内部联系在数量方面的反映，利用函数关系又可以对客观事物的规律性进行研究，因此如何寻求函数关系，在实践中具有重要意义。在许多问题中，不能直接找到所需的函数关系，但是根据问题所提供的情况，有时可以列出含有要找的函数及其导数的关系式，这样的关系式称为：微分方程。对其进行研究，找寻未知函数，称为解微分方程。本章主要介绍微分方程的一些基本概念和几种常用解法 微分方程的幂级数解法 当微分方程的解不能用初等函数或其积分式表达时，我们就要寻求其它解法。常用的有幂级数解法和数值解法。本节我们简单地介绍一下微分方程的幂级数解法。 求一阶微分方程（1）满 足初始条件的特解，其中函数 f (x , y)是、的多项式： . 这时我们可以设所特解可展开为 的幂级数 (2） 其中是待定的系数，把（2）代入（1）中，便得一恒等式，比较这恒等式 两端的同次幂的系数，就可定出常数 , 以这些常数为系数的级数（2）在其收敛区间内就是方程（1）满足初始条件 的特解。 例1求方程满足的特 解。 解这时，故设 ， 把及的幂级数展开式代入原方程，得 由此，比较恒等式两端x 的同次幂的系数，得 于是所求解的幂级数展开式的开始几项为 。 关于二阶齐次线性方程用幂级数求解的问题，我们先叙述一个定理： 定理如果方程（3）中的系数P(x)与Q(x)可在－R＜x＜R 内展开为x的幂级数那么 在－R＜x＜R内方程（3）必有形如 的解。 例 2 求微分方程的满足初始条件 , 的特解。 解这里在整个数轴上满足定理的条件。因此所求的解可在整个数轴上殿开成x的幂级 数（4） 由条件得。对级数（4）逐项求导，有 , 由条件得.于是我们所求方程的级数解及的形式已成为 （5） （6） 对级数（6）逐项求导，得 第四章 线性常微分方程的级数解法 4.1 常点邻域之级数解法 ① 常点邻域的级数解概念 ---- （二阶线性常微分方程的一般形式） 0)()(=+'+''w z q w z p w (4.1) ----（常点概念） 对于式(4.1)中,若)(z p 与 )(z q 在某点及其邻域内解析,则称此点为常点； 反之,若)(z p 与)(z q 至少一个在该点不解析,则称此点为奇点。 ----（常点邻域内解的存在定理） 若)(z p 与 ) (z q 在 R z z <-0内单值解析,则方程(4.1)在 R z z <-0内存在单值唯一的解析解。 ----（常点0z 邻域内之级数解的一般形式） 若 )(z p 与)(z q 在R z z <-0内单值解析,则对于式 (4.1),可设级数解∑∞ =-=0 0)(n n n z z a w ,再将 ) (z p 与 )(z q 在R z z <-0内展为泰勒级数,代入式(4.1)以 确定级数解之待定系数。 ② 勒让德方程之级数解 ----（勒让德方程形式） 0)1(2)1(2=++'-''-y l l y x y x (4.2) ----（在常点0=x 邻域内的级数解） 分析： 由1 2)(2-= x x x p 及2 1) 1()(x l l x q -+=，可知0=x 为常点；故可设：∑∞ ==0 n n n x a y ， 相应：∑∞ =-='1 1 n n n x na y ，∑∞ =--=''2 2)1(n n n x a n n y ， 代入方程(4.2)，得： )1(2)1()1)(2(0 2=++--- ++∑∑∑∑∞ =∞ =∞ =∞ =+n n n n n n n n n n n n x a l l x na x a n n x a n n ，即： n n a l l n n a n n )()1)(2(222--+=+++，或 n n a n n l n l n a ) 1)(2() 1)((2++++-=+；显然有： 02!2)1)((a l l a +-= ，13!3) 2)(1(a l l a +-=， 04! 4)12)(2)(1)((a l l l l a ++-+-=， 15! 5)4)(3)(2)(1(a l l l l a +-+-=，即 02)! 2() 12)(22()1)((a k l k l k l l a k +---+-= ， 012)! 12() 2)(12()2)(1(a k l k l k l l a k ++--+-= + ；相应级 数解为两个线性无关解的迭加： ∑∑∑∑∞ =++∞ =∞ =++∞ =+=+ = 1 21210 220 1 2120 22k k k k k k k k k k k k x A a x A a x a x a y (4.3)Mathematica函数大全(内置)
二阶线性常微分方程的幂级数解法
Mathematica的常用函数
Mathematica函数大全
高阶线性微分方程常用解法介绍
Mathematica常用指令
高阶方程的降阶法幂级数解法
二阶线性常微分方程的幂级数解法
Mathematica中的常用函数命令
Mathematical常用功能大全-精简版
高阶线性微分方程常用解法简介
微分方程的幂级数解法
线性常微分方程的级数解法