第6章 方程求根与解常微分方程
6.1实验目的
了解微分方程的通解、特解和近似解的概念。熟悉方程求根和常微分方程解的概念，熟悉Mathematica 软件的方程求根和求常微分方程解的命令，掌握用数学软件处理方程求根和常微分方程解的有关问题．
6.2实验准备
6.2.1数学概念
1．微分方程
2．微分方程的通解、特解
6.2.2数学软件命令
1. Solve[eqn, x]
功能：求多项式方程eqn 的所有根，当多项式方程的次数n ≤4时，给出eqn 所有根的准确形式, 当n>4时，不一定能求出所有的根, 此时，命令输出形式为
{ToRules[Roots[eqn, x ]]}
n 次多项式方程的一般形式为:
a 0+a 1x +a 2x 2+" +a n x n =0
式中a 0 ,a 1, a2,…,a n为常数。
2. Solve[{eqn1, eqn2, …, eqnk}, {x1, x2,…, xk}]
功能：求多项式方程组{eqn1, eqn2, …, eqnk}的所有根, 当其中每个多项式方程的次数n 4时, 给出所有根的准确形式, 否则，不一定能求出所有的根, 此时，命令输出形式为{ToRules[Roots[{eqn1, eqn2, …, eqnk}, {x1, x2,…, xk} ]]} 。
3. NSolve[eqn, x]
功能：求多项式方程eqn 的所有根的近似形式。
4. NSolve[{eqn1, eqn2, …, eqnk}, {x1, x2,…, xk}]
功能：求多项式方程组{eqn1, eqn2, …, eqnk}所有根的近似形式。
5. FindRoot[eqn, {x, x0}]
功能：求方程eqn 的在初值x0附近的一个近似根。
6. FindRoot[{eqn1,eqn2, ... }, {x, x0}, {y, y0}, ... ]
功能：求方程组{eqn1, eqn2, …}在初值(x0,y0,…)附近的一个近似根。
7.DSolve[eqn，y[x]，x]
功能：求常微分方程eqn 的通解y=f(x,C).
8.DSolve[{eqn，初始条件}，y[x]，x]
功能：求常微分方程eqn 满足初始条件的特解.
9.DSolve[{eqn1,eqn2, ... }, {y1[x],y2[x], ... }, x]
功能：求常微分方程组{eqn1,eqn2, ... }的解.
10. NDSolve[ eqns, y, {x, xmin, xmax}]
功能：求出自变量范围为[xmin, xmax]且满足给定常微分方程及初值条件eqns 的未知函数y 的数值解. 得到的解是以
{{未知函数名->InterpolatingFunction[range, ]}}
的形式给出的, 其中的InterpolatingFunction[range, ]是所求的插值函数表示的数值解, range就是所求数值解的自变量范围。
11.NDSolve[ eqns, {y1,y2, ... }, {x, xmin, xmax}]
功能：求出自变量范围为[xmin, xmax]且满足给定常微分方程及初值条件eqns 的未知函数y1,y2, ... 的数值解。
6.3实验任务
6.3.1基础实验
本实验熟悉数学软件命令操作。
1. 求下列代数方程的根
1）设 x 3 -4x 2 +9x - 10 = 0 ，求方程的所有根。
2）设x 2 -ax - 4b=0，求方程的所有根，a ，b 为常数。
3）求方程组
⎧x +3y =0 ⎨22⎩x +y =1
的所有根。
4）设 x 6 -x2 +2x - 3 = 0 ，求方程的所有根。
5）求方程组
⎧a 1+a 2=1⎪⎪x 1a 1+x 2a 2=14⎪⎪⎨21 2⎪x 1a 1+x 2a 2=9⎪⎪x 3a +x 3a =1
1122⎪16⎩
的所有根，这里x 1 ，x 2，a 1，a 2 是变量。
2．求下列超越方程的解
1）求方程组
⎧x =y 2
⎨⎩y =cos x
在(1,2)附近的根
2）求方程e x /2sin x =2在[0, 3]内的所有根。
3. 求解下列微分方程
1) 求常微分方程y ′=2ax 的通解，a 为常数
2) 求常微分方程y ′′+y =0的通解
3) 求常微分方程
⎧′x y ⎪y =+ y x ⎨⎪y (1)=2⎩
的特解。
⎧dx =y +1⎪⎧x (0)=−2⎪dt 4) 求方程组⎨且满足⎨的特解． (0)0dy y =⎩⎪=x +1⎪⎩dt
5) 求微分方程初值问y ′=x +y , y (0)=1在区间[0,1]内的数值解
6) 求y ′′−y ′=x ，满足y (0)=0，y ′(0)=0的特解．
6.3.2探索实验
本实验探索最速降线问题。
4. 一初始速度为零的质点，仅受到重力的作用，沿光滑固定的曲线由定点滑行到定点 (低于，但不在同一铅直线上) ．为使滑行的时间最短，确定该曲线具有什么形状的问题称为最速降线问题（如图6-1）。理论推导得出该最速降线可以由如下微分方程描述：
⎡1+y ′2⎦⎤=c 且y (0)=
0 y ⎣
图6-1
请给出该最速降线的曲线方程。
6.3.3应用实验
本实验研究传染病传播问题。
10小时后有2人被传染发病。5. 一艘游船载有1000人, 一名游客患了某种传染病，
由于这种传染病没有早期症状，故传染者不能被及时隔离。假设直升飞机将在50至60小时将疫苗运到，试估算疫苗运到时患此传染病的人数。
6.4实验过程
1．
1）In[1]: = Solve[x^3-4x^2+9x-10==0,x]
Out[1]= {{x -> 1 - 2 I}, {x -> 1 + 2 I}, {x -> 2}}
所以，所求全部根为 x1=1-2I，x2=1+2I ， x3=2， I 为虚数单位。
2）In[2]: = Solve[x^2-a*x-4*b==0,x]
⎧⎧11⎫⎧⎫⎫Out[2]= ⎨⎨x →a ⎬, ⎨x →a +⎬⎬ 22⎭⎩⎭⎭⎩⎩
所求全部根为
(
(x 1=11a −, x 2=a + 22(
(3）In[3]: = Solve[{x+3y==0,x^2+y^2==1},{x,y}]
⎧⎧⎧⎫Out[3]= ⎨⎨x →y →, ⎨x →y → ⎩⎭⎩⎩
4）In[4]: = Solve[x^6-x^2+2x-3==0 , x]
Out[4]= {ToRules[Roots[2 x – x2 + x 6 = = 3, x]]}
此结果说明求不出准确根, 改用NSolve 命令
In[5]: = NSolve[x^6-x^2+2x-3==0,x]
Out[5]= { {x -> -1.40825}, {x -> -0.465869 - 1.19413 I}, {x -> -0.465869 +
1.19413 I}, {x -> 0.608047 - 0.885411 I}, {x -> 0.608047 + 0.885411 I}, {x ->
1.12389} }
得所求全部6个近似根为
x1=-1.40825, x2=-0.465869- 1.19413 I, x3=-0.465869 + 1.19413 I,
x4= 0.608047 - 0.885411 I, x5= 0.608047 + 0.885411 I, x6= 1.12389 , I为虚数单位。
5）In[6]:=Solve[{a1+a2==1,
x1*a1+x2*a2==1/4,
x1^2*a1+x2^2*a2==1/9,
x1^3*a1+x2^3*a2==1/16},
{a1,a2,x1,x2}]
11Out[6]=⎧212−,a2→212+⎨a1→424424⎩(
(
x2→11⎫15−, x1→15+⎬ 4242⎭(
(In[7]: = N[%] （*显示根的近似形式*）
Out[7]={{a1 -> 0.281461, a2 -> 0.718539, x2 -> 0.112009, x1 -> 0.602277}, {a1 -> 0.718539, a2 -> 0.281461, x2 -> 0.602277, x1 -> 0.112009}
In[8]: = N[%, 8] （*显示根取8位有效数字的近似值*） Out[8]= {{a1 -> 0.28146068, a2 -> 0.71853932, x2 -> 0.11200881, x1 ->
0.60227691}, {a1 -> 0.71853932, a2 -> 0.28146068, x2 -> 0.60227691, x1 -> 0.11200881}}
2.
1）In[1]: = FindRoot[{x==y^2,y==Cos[x]},{x,1},{y,2}]
Out[1]= {x -> 0.641714, y -> 0.801071}
所求的根为
x = 0.641714, y = 0.801071。
2) In[2]: = Plot[Sin[x]*Exp[x/2]-2,{x,0,3}]
（*画出方程对应的函数在[0,3]内的图形*）
输出图形为
Out[2]= -Graphics- （*图中可知方程在1.5和2.5附近有根*） In[3]: = FindRoot[Sin[x]*Exp[x/2]-2==0,{x,1.5}]
（*求方程在1.5附近的根*）
Out[3]= {x -> 1.41171}
In[4]: = FindRoot[Sin[x]*Exp[x/2]-2==0,{x,2.5}]
（*求方程在2.5附近的根*）
Out[4]= {x -> 2.54786}
所求根为x1 = 1.41171 , x2= 2.54786。
3.
1)In[1]:= DSolve[y'[x]==2a*x,y[x],x]
Out[1]= {{y[x] -> a x 2 + C[1]}}
即得本问题的通解 y =c 1+ax 2
2) In[2]:= DSolve[y''[x]+y[x]==0,y[x],x]
Out[2]= {{y[x] -> C[2] Cos[x] - C[1] Sin[x]}}
即得本问题的通解: y =c 2cos x −c 1sin x
C 1，C 2是任意常数。
3) In[3]:= DSolve[{y'[x]==x/y[x]+y[x]/x,y[1]==2},y[x],x]
Out[3]= {{y[x] -> Sqrt[x2 (4 + 2 Log[x])]}}
本问题的特解为:
y =
下面的操作可以检验所求特解的正确性
In[4]:= y[x_]:= Sqrt[x^2*(4+2Log[x])] （*将所求特解定义为一个函数*） In[5]= y[1] （*检验初始条件y[1]=2*）
Out[5]=2 （*初始条件成立*）
In[6]:= D[y[x],x]-x/y[x]-y[x]/x; （*将微分方程左端项与右端项相减*） In[7]:= Simplify[%] （*化简计算结果*）
Out[7]= 0 （*说明所求特解正确*）
4）In[8]:= DSolve[{x'[t]==y[t]+1,
y'[t]==x[t]+1,
x[0]== -2,
y[0]==0},
{x[t],y[t]},t]
Out[8]:=x [t ]→−e −t (1+e t ), y [t ]→−e −t (−1+e t )
5) In[9]:= NDSolve[{y'[x]==x+y[x],y[0]==1},y,{x,0,1}]
Out[9]= {{y -> InterpolatingFunction[{0., 1}, ]}}
上面显示的是所求数值解的替换形式，为得到本问题数值解，再键入: {{}}In[10]:= f=y/.%[[1]]
（*把InterpolatingFunction[{0., 0.5}, ]赋值给变量f*）
Out[10]= InterpolatingFunction[{0., 1}, ]
In[11]:= Plot[f[x],{x,0,1}] （*画出所求未知函数数值解的图形*）
输出图形为
Out[11]= -Graphics-
In[12]:= Table[f[x],{x,0,0.5,0.1}]
（*计算数值解在0,0.1,0.2,0.3,0.4,0.5的值*）
Out[12]= {1., 1.11034, 1.24281, 1.39972, 1.58365, 1.79744}
6) In[13]:= DSolve[{y''[x]-y'[x]==x,y[0]==0,y'[0]==0},y[x],x] 1⎧⎧⎫⎫ Out[13]=⎨⎨y [x ]→(−2+2e x −2x −x 2)⎬⎬ 2⎭⎭⎩⎩
4.
1）分析
2′⎤1+y 由于所给微分方程y ⎡⎣⎦=c 不能直接用DSolve 解出，借助换元法将其转化
dy 1⎡1+y ′2⎦
⎤=c 可解出为参数方程的形式。令y =c (1−cos 2t )，再由y ⎣= 2dx dy
dx dx dy 于是，然后积分得到x =x (t ) ，即得到一个参数解．程序如下： ==dt dy dt dy
dx
2）实验操作
In[1]:= y=c*(1-Cos[2*t])/2；
yt=D[y,t]；
yx=Sqrt[c/y-1]；
xt=yt/yx；
x=Integrate[xt,t]；
Print["x=",x]
Print["y=",y]
输出的结果为：
c −Cos[t ]2+t Cot[t ]
1 y=c (1−Cos[2t ]) 2
代入y (0)=0，可以得到c 1=0．令a =c ，θ=2t ，则有最速降线方程为 2
⎧⎪x =a (θ−sin θ) ⎨⎪⎩y =a (1−cos θ)
其中θ为参量，a 为待定参数。从方程可以知道这是高等数学中的摆线标准参数方程．对于具体的最速降线，只要给出B 点的坐标，就可以确定出c 的值，从而最速降线也就完全确定了。
5．
1) 问题分析
设y(t)为发现第一个病人后t 小时时刻的传染人数，则有y(t)对时间t 的导数dy 可以描述该传染病的传染速率。常识表明传染病的传染速率既受到传染人数dt
的影响又受未被传染人数的影响。一般情况下有染病人数越多传染速度越快（因为有很多的传染源）、未被传染人数越多传染速度也越快（因为会有很多的人染病）。因此其影响关系都为正比关系，本题中在t 时刻未被传染的人数为1000-y(t),于是我们可以用如下微分方程描述传染速率：
dy =ky (1000−y ), dt y (0)=1, y (10)=2，其中k 为比例常数
求解此微分方程即可。
2) 实验操作
In[1]:=DSolve[y'[t]==k*y[t]*(1000-y[t]),y[t],t] ⎧⎧⎫1000e 1000kt ⎫⎪⎪ Out[1]:= ⎨⎨y [t ]→1000kt 1000C [1]⎬⎬e −e ⎪⎭⎪⎩⎩⎭
1000e 1000kt
因此得到解y [t ]=1000kt 。 1000C [1]e −e
由y (0)=1可得c =999, 由y (10)=2，可以得1000k=0.0694149
In[2]:= y1[t_]:=1000/(1+999*Exp[-0.0694149t])
于是得到t 小时时刻的传染人数
y (t ) =1000 −0.0694149t 1+999e
In[3]:= Plot[y1[t],{t,0,200}]
输出图形为
Out[3]:=-Graphics-
In[4]:= y1[50]
Out[4]:= 31.1888≈32
In[5]:= y1[60]
Out[5]:= 60.548≈61
因此在t=50小时患此传染病的人数约为32人，在t=60小时患此传染病的人数约为61人。
从这些数字可以看到从50小时到60小时这10小时之间被传染发病的人数几乎翻了一倍，因此在传染病流行期间应该及时采取措施是很重要的。如果不采取措施，通过y(t)的图形可以看到当t 在50到150小时之间传染最快，且当t 趋于无穷大时，y(t)是趋于1000人，即导致全游艇人数都被传染。
6.5思考与提高
1. Solve命令和FindRoot 命令有什么区别？
2. 如果超越方程中含有自变量之外的字母，还能用FindRoot 命令求根吗？ 3. 求常微分方程的DSolve[eqn, y[x]，x] 命令中的“y[x]，x”能否写为“y[t]，t”？
6.6练习内容
1. 求下列代数方程的根
1）设 2x3 +x 2 - 3 = 0 ，求方程的所有根。
2）设abx 2+(a4+b4)x +a3b 3=0，ab ≠0, 求方程的所有根，a ，b 为常数。
3）求方程组
⎧2x −3y =5 ⎨22⎩4x −9y =15
的所有根。
4）设 x 8 –x4 +2x - 3 = 0 ，求方程的所有根。
2. 求下列超越方程的解
1) 求方程x sin x – 0.5=0在[-5, 5]内的所有根.
2) 求方程cos3x+2cosx=0在[-10, 10]内的所有根。
3. 求解下列微分方程
1) 求y ′cos 2x +y =tan x 的通解．
2) 求y ′+y +xy 2=0的通解．
3) 求常微分方程组y ′(t ) =x (t ) ，x ′(t ) =y (t ) 的通解
dy 2−3x 2
y =1，且满足y(1)=0的特解． 4) 求+3dx x
⎧x ′+5x +y =e t ⎧x (0)=1，满足⎨的特解． 5) 求⎨⎩y (0)=0⎩y ′−x −3y =0
4. 设有一个弹性系数为的弹簧竖直悬挂着，它的上端固定，下端系上质量为的物体。如果在初始时刻物体运动到最低点，运动距离为，求物体的运动方程。 5. 一敌舰在某海域内沿着正北方向航行时，我方战舰恰好位于敌舰的正西方向1n mile处．我舰向敌舰发射制导鱼雷，敌舰速度位0.42n mile/min，鱼雷速度位敌舰速度的2倍。试问敌舰航行多远时将被击中？
6. 选择一道涉及常微分方程求解的应用题做之，对解答过程中出现的需要求解的数学式子用Mathematia 数学软件命令来求解。