微分方程[1].ppt

微分方程建模微分方程建模适应范围适应范围 与变化率有关的各种实际问题应用三步曲应用三步曲 (1)建模 即根据实际问题建立起适当的微分方程，给出其定解条件.(2)求解 求出所建立的微分方程的解 (3)翻译 用所得结果来解释一些现象，或对问题的 解决提出建议或方法建议:模型要详略得当 在用微分方程解决实际问题的过程中一定要意识到实际问题是十分复杂的，微分方程只能是在一定程度上对问题的一种近似描述，只要结果的误差在一定范围内即可.任何模型都不可能把影响问题的所有因素都反映在微分方程中，或者要求所得结果十分精确.一个好的微分方程模型是在实际问题的精确性和数学处理的可能性之间的一个平衡.线性微分方程的应用举例线性微分方程的应用举例在电感上的电压降为 由Kirchhoff回路电压定律知：沿着任一闭合回路的电压降的代数和为零。我们得到电流 所满足的微分方程为：例例1 1：RL串联电路由电阻、电感、关闭合后电路中的电流强度 电源组成的串联电路，求开解解：当电路中电流为 时，在R上的电压降为 取开关闭合时刻为0，则 又 将初始条件 代入齐次方程通解得 故当开关闭合后，电路中的电流强度为：求得齐次方程通解为是方程特解，例例2 2 湖泊的污染湖泊的污染 设一个化工厂每立方米的废水中含有3.08kg盐酸，这些废水流入一个湖泊中，废水流入的速率20立方米每小时.开始湖中有水400000立方米.河水中流入不含盐酸的水是1000立方米每小时,湖泊中混合均匀的水的流出的速率是1000立方米每小时,求该厂排污1年时,湖泊水中盐酸的含量。解:设t时刻湖泊中所含盐酸的数量为考虑 内湖泊中盐酸的变化。例例2 2 湖泊的污染湖泊的污染设一个化工厂每立方米的废水中含有3.08kg盐酸，这些废水流入一个湖泊中，废水流入的速率20立方米每小时.开始湖中有水400000立方米.河水中流入不含盐酸的水是1000立方米每小时,湖泊中混合均匀的水的流出的速率是1000立方米每小时,求该厂排污1年时,湖泊水中盐酸的含量。解解:设t时刻湖泊中所含盐酸的数量为考虑 内湖泊中盐酸的变化。因此有该方程有积分因子两边同乘以后,整理得积分得利用初始条件得初值问题初值问题 的解为的解为 初值问题初值问题 的解为的解为 雪球融化问题雪球融化问题设雪球在融化时体积的变化率与表面积成比例，且融化过程中它始终为球体，该雪球在开始时的半径为6cm，经过2小时后，其半径缩小为3cm。求雪球的体积随时间变化的关系。解:设t时刻雪球的体积为，表面积为，球体与表面积的关系为 由题得引入新常数 再利用题中的条件得分离变量积分得方程得通解为再利用条件 确定出常数C和r代入关系式得 t的取值在 之间。有一段时间,美国原子能委员会(现为核管理委员会)是这样处理浓缩放射性废物的,他们把这些废物装入密封性能很好的圆桶中,然后扔到水深300英尺的海里。这种做法是否会造成放射性污染,很自然地引起了生态学家及社会各界的关注。原子能委员会一再保证,圆桶非常坚固,决不会破漏,这种做法是绝对安全的。然而一些工程师们却对此表示怀疑,他们认为圆桶在和海底相撞时有可能发生破裂。而原子能委员会有专家们则仍然坚持自己的看法。于是,双方展开了一场笔墨官司。究竟谁的意见正确呢?看来只能让事实说话了。问题的关键在于圆桶到底能承受多大速度的碰撞,圆桶和海底碰撞时的速度有多大?放射性废物的处理问题放射性废物的处理问题 大量破坏性实验,发现圆桶在40英尺秒的冲撞下会发生破裂,剩下的问题就是计算圆桶沉入300英尺深的海底时,其末速度究竟有多大了。美国原子能委员会使用的是55加仑的圆桶,装满放射性废物时的圆桶重量为W527.436磅,而在海水中受到的浮力B470.327磅。此外,下沉时圆桶还要受到海水的阻力,阻力Dv,其中C为常数。工程师们做了大量实验,测得C0.08。现在,取一个垂直向下的坐标,并以海平面为坐标原点(0)。于是,根据牛顿第二定律建立圆桶下沉时应满足方程 质量质量加速度加速度=重力重力-浮力浮力-摩擦阻力摩擦阻力 模型及其解oymgBD困难：无法知道下沉到海底的时间积分和代入初始条件得：最后再用数值计算可以得到水深300时的速度大小。借助数值方法求出v(300)的近似值。计算结果表明,v(300)45.1英尺秒40英尺秒。工程师们的猜测是正确的,他们打赢了这场官司。现在,美国原子能委员会已改变了他们处理放射性废物的方法,并明确规定禁止将放射性废物抛入海中。：有高为1米的半球形容器,水从它的底部小孔流出,小孔横截面积为1平方厘米(如图).开始时容器内盛满了水,求水从小孔流出过程中容器里水面的高度h(水面与孔口中心间的距离)随时间t的变化规律.例例 解：解：由力学知识得,水从孔口流出的流量为流量系数孔口截面面积重力加速度小孔流水问题设在微小的时间间隔水面的高度由h 降至 ,比较(1)和(2)得:即为未知函数的微分方程.可分离变量所求规律为值得进一步探讨的问题:漏斗型的容器 由于水的张力的原因，每次水都无法全部留尽，总会剩一小部分在容器中。如何才能让水尽可能少的留在容器中？我们知道，水与容器接触的面积越大，留在容器中的水就越多先讨论一下漏斗的模型。y 容器的位置容器的位置 可否将容器倾斜，使上部的面积大于下部的面积，使水流的速度更快？倾斜角度？容器的运动状态容器的运动状态 容器的运动状态对流水的速度是肯定会造成影响的，考虑极限的状态，如果容器以大于等于当地重力加速度的加速度竖直向下运动，那么，容器里的水就不会流出。容器以不同的方式运动时对水的流出时间有多少影响？有没有一种运动状态能加快水流的速度呢？涡流的影响涡流的影响 涡流对水流的速度是有一定影响的。拿一个水桶反复做这样的试验：首先将桶装满水，记录水面的高度，然后拔出塞住孔口的塞子，让水自然从桶破了的孔中流出，测量流出的时间，然后反复从同一高度作相同的试验，最后求出水自然流尽所需时间的平均值；然后从同一高度作相同的试验，不同的是用一根棍子绕同一方向在水中搅动，使其产生涡流，然后重复上面的步骤。最后发现通过两种方法测得的水流尽所需时间的平均值有较大的差距，于是猜想有无涡流或许对水流的速度也是有一定影响的。疾病的传播疾病的传播 SARS的预测 SARS是21世纪第一个在世界范围内传播的传染病。SARS从2002年11月份开始在我国和世界范围内流行，到2003年6月23日为止，世界卫生组织(WHO)报道的SARS患者已经达到了8459人，其中802人死亡。中国是SARS流行的重灾区，到2003年6月23日为止的SARS患者为5326人，其中347人死亡。给人民生活和国民经济的发展带来了巨大的影响。SARS是由一种冠状病毒引起的传染性很强的呼吸道传染病，它主要通过近距离空气飞沫以及接触病人呼吸道分泌物和密切接触进行传播，也可能通过病人飞沫污染物、如通过手、衣物、食物、水或环境等途径传播。SARS潜伏期一般为2-11天，在潜伏期无感染。SARS患者的主要症状有：发热（体温38以上）为首发症状，多为高热，并可持续1-2周以上，可伴有寒战或其他症状，包括头痛、全身酸痛和不适、乏力，部分病人在早期也会有轻度的呼吸道症状(如咳嗽、咽痛等）。SARS患者治愈后不会再被感染。假设：假设：1）单位时间感染的人数与现有的感染者成比例；单位时间感染的人数与现有的感染者成比例；2）单位时间内治愈的人数与现有的感染者成比单位时间内治愈的人数与现有的感染者成比 例；例；3）单位时间内死亡的感染者人数与现有的感染者 成比例；4）SARS患者治愈恢复后不再被感染；5）各类人口的自然死亡可以忽略；6）忽略迁移的影响。令I(t)是第t天时SARS感染者的数量，b(t)为感染率，d(t)为死亡率，c(t)为治愈率模型根据当时SARS数据是按天公布这个特点,离散化得 I(t+1)=I(t)+r(t)I(t)只要知道开始时SARS的感染人数和函数r(t)，就可以利用该模型进行预测。r(t)的估计 r(t)=I(t+1)-I(t)/I(t)利用实际数据计算，再进行曲线拟合服药问题服药问题 医生给病人开处方时必须注明两点:服药的剂量和服药的时间间隔.超剂量的药品会对身体产生严重不良后果,甚至死亡,而剂量不足,则不能达到治病的目的.已知患者服药后,随时间推移,药品在体内逐渐被吸收,发生生化反应,也就是体内药品的浓度逐渐降低.药品浓度降低的速率与体内当时药品的浓度成正比.当服药量为A、服药间隔为T,试分析体内药的浓度随时间的变化规律.体内药的浓度随时间的变化规律井深的估算井深的估算假如你站在井边且身上带着一只具有秒表功假如你站在井边且身上带着一只具有秒表功 能的计算器，你也许会出于好奇心想用扔下能的计算器，你也许会出于好奇心想用扔下 一块石头听回声的方法来估计井的深度，一块石头听回声的方法来估计井的深度，假定你能准确地测定时间，你又怎样来推算假定你能准确地测定时间，你又怎样来推算 井的深度呢，请你分析一下这一问题。井的深度呢，请你分析一下这一问题。我有一只具有秒我有一只具有秒 表功能的计算器。表功能的计算器。方法一方法一假定空气阻力不计，可以直接利用自由落体运动的公式假定空气阻力不计，可以直接利用自由落体运动的公式来计算。例如，来计算。例如，设设t=4=4秒，秒，g=9.81=9.81米米/秒秒2 2，则可求得，则可求得h78.578.5米。米。我学过微积分，我可以做我学过微积分，我可以做 得更好，呵呵。得更好，呵呵。除去地球吸引力外，对石块下落影响最大的当属空气阻力。根据流体力学知识，此时可设空气阻力正比于石块下落的速度，阻力系数K为常数，因而，由牛顿第二定律可得：令k=K/m，解得 代入初始条件 v(0)=0，得c=g/k，故有 再积分一次，得：若设k=0.05并仍设 t=4秒，则可求 得h73.6米。听到回声再按跑表，计算得到的时间中包含了 反应时间 进一步深入考虑进一步深入考虑不妨设平均反应时间 为0.1秒，假如仍 设t=4秒，扣除反应时间后应 为3.9秒，代入 式，求得h69.9米。多测几次，取平均值再一步深入考虑再一步深入考虑代入初始条 件h(0)=0，得到计算水井高度的公式：将e-kt用泰勒公式展开并 令k 0+，即可得出前面不考虑空气阻力时的结果。还应考虑回声传回来所需要的时间。为此，令石块下落 的真正时间 为t1，声音传回来的时间记 为t2，还得解一个方程组：这一方程组是非线性的，求解不太容易，为了估算井深竟要去解一个非线性主程组似乎不合情理 相对于石块速度，声音速度要快得多，我们可 用方法二先求一次 h，令t2=h/340，校正t，求石块下落时间 t1t-t2将t1代入式再算一次，得出井深的近似值。例如，若h=69.9米，则 t20.21秒，故 t13.69秒，求得 h62.3米。速度V运动，方向永远指向P点,求M点的运动在轴上有一点P以常速度a沿着轴追线问题追线问题平面上另有一点M,它以常正向移动;在轨迹.解:首先我们建立点首先我们建立点M M运动时所满足的微分方程运动时所满足的微分方程模型模型标,根据条件有 图3.1标,根据条件有（3.1.7）以记点M在时刻t的坐标，以X记图3.1点P在时刻t的横坐标，表示P点在t=0的横坐（3.1.6）把（3.1.6）代入（3.1.8），并记（3.1.8）得:上式两边关于作为自变量，把求导得（3.1.9）即由（3.1.9）和（3.1.10）得到M的追线方程 于是（3.1.11）变为 又由得:（3.1.10）（3.1.11）令则讨论由图3.1知,在点M未追上点P之前,点P的横坐标总大于点M 的横坐标,即当时所以积分上式得:即由此得（3.1.12）从 得解,即点M沿Ox轴移动.假设开始追逐是,点 P 和 M 同在一条平行于轴的直线上,并记它们的位置为及显然有,由此得,从而得:由上式得:（3.1.13）（3.1.14）(3.1.13)减去(3.1.14)得:即为了使点M有可能追上点P,我们假设（3.1.15）此时,由(3.1.15)得到追线方程为:由于初始点在追线上,即当时,因此得:从而得追线方程:当时,就得到相遇点的坐标是追上所需的时间是悬链线问题悬链线问题有一绳索悬挂在A和B两点(不一定是在同一水平线),如图3.2所示.设绳索是均匀的,柔软的,仅受绳本身的重量作用,它弯曲如图中的形状,试确定该绳索在平衡状态时的形状.解解:设C是其最低点,选取坐标系如图中所示,且轴通过C点.A AB BC CO O图图3.23.2考虑绳索在最低点C与点之间的一段,这一段在下面三个力的作用下平衡:(1)在点P的张力T,方向沿着P点的切线方向;(2)在点C的水平张力H;(3)CP段的垂直的重量,记为,设它作用在某一点Q处,不一定是CP的中心,见图3.3由于平衡关系,这些力在轴(水平)方向的代数和为0,在轴(垂直)方向的代数和也必须为0.T TQ QC CH H图图3.33.3现将张力分解为两个分力：水平方向分，垂直方向分力为此时，在力为轴方向向左而向右；在轴方向，向下而向上，按平衡关系有：两式相除，并利用关系式（为点的切线斜率）得：是在最低点处的张力，是常数，但依赖于，将上式两边对微分得其中表示在水平方向上，每增加单位距离时，段弧所增加的重量设绳索的密度为，则有其中表示从点算起的弧长，我们需要求出，因为（3.1.1）或又由于故从而方程（3.1.16）化为：（3.1.1）记为绳索最低点到坐标原点的距离，则有：（3.1.1）是一个不显含自变量的方程，令则方程（3.1.1）化为分离变量，积分得：即式中把初始条件代入（3.1.1）上式得：（3.1.1）上式两端同时乘以，故（3.1.1）变为得：两式相加得：即积分上式，得：把初始条件代入上式得为简单起见，假设此时，从而得绳索的方程：上式表示的曲线叫做悬链线此时，绳索在最低点与点之间的一段弧长是：Volterra 捕食捕食-被捕食模型被捕食模型 设有捕食种群和食饵种群生活在同一小环境中,建立微分方程组来研究两种群个体数量随时间的变化趋势.设 t 时刻食饵和捕食者的数量或密度分别为假设个体不区分大小,而且没有个体向环境输入或从环境输出,当环境中不存在捕食者时,食饵种群的增长规律用下述Logistic方程来描述上式左端表示被捕食者的相对增长率;右端的常数称为内禀增长率,为环境的容纳量,由(4.1.6)可以看出,(4.1.6)因此当时,种群规模增长,时,种群规模减小.反映了环境能保证食饵个体数量变化时最合适的容量,把(4.1.6)改写形式是其出生率 减去死亡率(4.1.7)其中项反映了以下事实:即在容纳量一定的条件下,的增大,将使每一个体平均的生活条件降低,从而影响种群的相对增长率,因此 或称为密度制约项.由于捕食者的存在,将使食饵的增长率减少,设单位总量成正比,注意到 t 时刻有y(t)个捕食者,它们在时间内每个捕食者吃掉的食饵数量与该时刻食饵的单位时间内吃掉食饵的总数量应为为常数,对于捕食种群,当不存在食饵种群时,仍用Logistic于是(4.1.7)变为方程来描述增长规律,即当存在食饵种群时,被捕食者吃掉的食饵将转化为能量去生育后代,设转化系数为则捕食种群的增长规律为其中式中项反映了捕食者仅以食饵 为生.这样我们得到一个Volterra 捕食-食饵系统(4.1.10)池水含盐量问题的推广 现在有3个水池，盛有不含盐的水各 1000 ，从t=0时刻起以每分钟4 /分的速率向第一个水池内注入盐水，同时又以每分钟3 /分的速率从该水池中流出搅拌均匀的盐水。第一个水池中流出的盐水进入第2个水池，同时再以2(1-cost)/分的速率向第2个水池中加入盐水，第2个水池中搅拌均匀的盐水又以2 /分的速率流入第3个水池，同时再以2(1-sint)/分的速率向第3个水池中加入盐水，第3个水池中搅拌均匀的盐水又以1 /分的速率出。已知输入盐水的含盐率为每立方米1千克，用计算机仿真这3个水池里盐水的变化过程，计算各水池中水的体积、含盐量、含盐率的变化。池水含盐量推广问题的示意图池水含盐量推广问题的示意图池水含盐量推广问题的微分方程模型池水含盐量推广问题的微分方程模型 令 分别是3个池子中t时刻的含盐量,分别是3个池子中t时刻水的体积,则restart;sys1:=diff(w1(t),t)=1,diff(w2(t),t)=1+2*(1-cos(t*Pi/180),diff(w3(t),t)=1+2*(1-sin(t*Pi/180);dsolve(sys1,w1(0)=1000,w2(0)=1000,w3(0)=1000,w1(t),w2(t),w3(t);assign(%);w1(t);plot(w1(t),w2(t),w3(t),t=0.200);dsolve(diff(x1(t),t)=16-x1(t)/w1(t),x1(0)=0,x1(t);assign(%),x1(t);dsolve(diff(x2(t),t)=3*x1(t)/w1(t)+8*(1-cos(t*Pi/180)-2*x2(t)/w2(t),x2(0)=0,x2(t);assign(%),x2(t);dsolve(diff(x3(t),t)=2*x2(t)/w2(t)+8*(1-sin(t*Pi/180)-x3(t)/w3(t),x3(0)=0,x3(t);assign(%);x3(t);plot(x1(t),t=0.200);水的变化水的变化池1-红，池2-黄，池3-绿池水含盐量推广问题的计算流程池水含盐量推广问题的计算流程 推广问题比较复杂，理论分析很困难，数值计算时难度相当，只是计算量有所增加输入初始值，取定步长在每一个循环内 计算各池中水的体积 含盐量和含盐率输出需要的结果数值计算程序数值计算程序:Euler:Euler 折线法折线法(Maple)(Maple)restart:h:=1.;t1:=h:s11:=0.:w11:=1000.:r11:=s11/w11:s21:=0.:w21:=1000.:r21:=s21/w21:s31:=0.:w31:=1000.:r31:=s31/w31:for i from 2 to 200 do ti:=i*h;s1i:=s1i-1+4*h-3*h*r1i-1;w1i:=w1i-1+h;r1i:=s1i/w1i;s2i:=s2i-1+r1i-1*3*h+2*(1-cos(ti*3.14/180)*h-2*h*r2i-1;w2i:=w2i-1+h+2*(1-cos(ti*3.14/180)*h;r2i:=s2i/w2i;s3i:=s3i-1+r2i-1*2*h+2*(1-sin(ti*3.14/180)*h-h*r3i-1;w3i:=w3i-1+h+2*(1-sin(ti*3.14/180)*h;r3i:=s3i/w3i;end do:p1:=seq(tn,s1n,n=1.200):p2:=seq(tn,s2n,n=1.200):p3:=seq(tn,s3n,n=1.200):q1:=seq(tn,r1n,n=1.200):q2:=seq(tn,r2n,n=1.200):q3:=seq(tn,r3n,n=1.200):plot(p1,p2,p3);plot(q1,q2,q3);含盐率的变化含盐率的变化池1-红，池2-绿，池3-黄微分方程组应用举例微分方程组应用举例 微分方程组在工程技术中的应用是非常广泛的，涉及的领域有机械、电工技术、通讯、医学等许多领域。下面给出几个例子，用以阐述关于微分方程组的实际应用及其建模思想。1 1、两自由度的振动问题两自由度的振动问题 常系数线性方程组在工程技术与科学研究中有很多应用，不少问题都归结为它的求解问题。例如两个自由度的振动问题。为了消除不需例如两个自由度的振动问题。为了消除不需要的振动，常常在振动系统中设置减振器。要的振动，常常在振动系统中设置减振器。其中的受力情况。是原机械部件的质量；是减振器的质量；和是两个弹簧，它们的弹性系数（或称为刚度）也分别用和是减速器表示；（假定阻力与速度成正比）的阻尼系数；是强迫力；和分别表示和距它们的平衡位置的位移。体和下面来建立这个系统的运动方程。先分别考虑物律，有方程(1)物体的受力情况假定弹簧和都满Hook定律。当物体位移时，物体同时位移同时位移这时，弹簧变形（拉长或压缩）的长度为。因此，这时弹簧的弹性力是（力的方向与位移方向相反），由牛顿第二定(2)物体 的受力情况1)沿位移方向的外力：2)阻尼力：（方向与速度方向相反）；3)这时物体 受到两个弹簧的作用：弹簧 的弹性力 ；弹簧 的弹性力 。由牛顿第二定律，有因此，上述运动系统满足微分方程组上述方程组在变换 及 之下，就变成了一个常系数线性非齐次方程组 (4.6.1)为了简单，只考虑无阻尼自由振动的情形，即 。于是，(4.6.1)变成 (4.6.2)它的特征方程为 设 ，则(4.6.3)可写成解之，得容易看出，因此，可令 ，这时因此，方程(4.6.3)最后可写成故(4.6.3)的特征根全是纯虚根。(4.6.1)解的第一个函数可写成形如部件的运动频率。这个事实可以用来防止机器或这是两个简谐运动的叠合。每一个简谐运动的角频（即与）均与减振器的参数与有关。因此，调整与可以改变原设备与外力发生共振现象，以及减轻外力的干扰等。在研究传染病的传播、种群生态、环境的污染、药物在人体内的分布等问题中，经常把所研究的事物看成由有限个部分组成的系统，而每个部分称为一个仓室仓室(compartment)。它具有以下特点：(1)每个仓室有固定的容量，内含每个时刻都均匀分布着的物质（或能量）。(2)各个仓室间以及仓室与外部环境间均可进行物质（或能量）交换，并服从物质（或能量）守恒定律。这样的系统称为仓室系统仓室系统。下面根据人体内胆固醇的吸收、合成、排泄等机理来建立胆固醇流动仓室模型。2 2、胆固醇流动的仓室模型、胆固醇流动的仓室模型 人体内血浆中的胆固醇既可从食物中吸取，也可以从体内通过肝脏等器官合成，且血浆中过量的胆固醇通过肝脏、肠等排出体外。把具有吸收和排泄胆固醇功能的器官组成的系统称仓室 ，如血浆、肝脏、肠道和皮肤等都属于仓室 。把与胆固醇有关的其他器官构成的系统称为仓室 ，如动脉壁等属于仓室 ，即在仓室 中，胆固醇既不能从外部吸收，也不能向外部排放，两个仓室都可以合成胆固醇，且两个仓室的胆固醇可以相互交换如图4.4所示图图4.44.4图 4.4 给出了两个仓室中胆固醇的转移率及 ，其中 代表单位时间每克胆固醇从仓室 流动到仓室 的量，代表单位时间每克胆固醇从仓室 流动到仓室 的量。其他参数意义为：代表单位时间内每克胆固醇从仓室 排除的量。代表单位时间通过从食物中吸收而输进仓室 的胆固醇的量。是代表单位时间通过人体合成而流入仓室 的胆固的量。(4.6.4)设和分别代表在时刻，仓室和中胆固醇的总量。通过上述分析和假设，可建立下面关于胆固醇的仓室模型为求解方程组(4.6.4)，将(4.6.4)写成向量形式 ，即这里 。方程组(4.6.4)对应的齐次方程组的特征方程为其判别式为因此 ，故系数矩阵 有两个负实部特征根的通解为其相应的特征向量为因此，是方程的基解矩阵，又因为是方程的解。故方程直接计算 得这样可得方程组(4.6.4)的通解为3 3、人造卫星的轨迹方程、人造卫星的轨迹方程我们知道，人造卫星在最后一段时间运载火箭熄灭之后，即进入它的轨道，轨道的形状因发射角度和发射速度的不同，而分别出现椭圆、抛物线或双曲线。下面来讨论这些问题。地球与人造卫星是相互吸引的，但因二者的质量相差很大，因此，可假设地球是不动的，又因人造卫星的体积与地球相比是很小的，故可把它看作质点。为简单起见，我们不考虑太阳、月亮和其它星球的作用，并略去空气阻力。在上面的假设下，可把问题归结为：从地球表面上一点 ，以倾角 ，初速度 射出一质量为 的物体，如图4.5所示，求此物体运动的轨道方程。图4.5 根据万有引力定律，地球对卫星的引力大小为过发射点和地心的直线作轴，轴与发射方向所成的平面为平面，平面通过地心，取垂直于轴且过地心的直线为轴，取开始发射时间为，经过时间后，卫星位于点，下面建立和所满足的方程。其方向指向地心，其中 是引力系数，是地球质量，是地球与卫星间的距离。如图4.6所示，这个引力在 轴方向上的分力分别为卫星在 轴上所获得的分加速度分别为 和 。由牛顿第二定律，得到卫星的运动方程为：轴上的分量分别为当时，卫星在地表面以倾角，初速度射出，所以，在时，。卫星的初速度在因此，初始条件为 （4.6.6）下面利用首次积分法来求方程组(4.6.5)满足初始条件(4.6.6)的解。将方程(4.6.5)两端消去 后，以 乘第一个方程，以 乘第二个方程。然后相减，得因为故有对上式两边积分，得首次积分为再以 乘方程组(4.6.5)中第一个方程，以 乘第二个方程，然后两式相加，得由于及从而得积分得上式另一首次积分于是，原方程组(4.6.5)降为较低阶的方程组 将它们代入(4.6.7)得作极坐标变换，并求导 两式联立消去 ，得 这就是卫星运动轨道的极坐标参数方程。若将参这里得到一个仅含有一个未知数得一阶微分方程，若由此解出，代入(4.6.9)中第一式，便可确定，由此得到数消去，便得出卫星运动轨道的极坐标方程。利用分离变量求该方程的解得为此，由(4.6.9)的第一式求得,并代入(4.6.10)得，则上式化为整理得令知或(4.6.11)这就是所求的卫星运动轨道的极坐标方程，其中有三个任意常数(或)，它们可由初始条件(4.6.6)确定。注意到当时，因此，。且由(4.6.6)及(4.6.8)把它们代入(4.6.9)及(4.6.11)得 (4.6.12)我们知道，(4.6.11)式是圆锥曲线的极坐标方程当 时，轨道是圆；当 时，轨道是椭圆；当 时，轨道是抛物线；当 时，轨道是双曲线 下面进一步讨论卫星发射的初速度与卫星轨道形状的关系因为故上故将(4.6.12)式中的代入上式，并整理得(4.6.13)注意到(4.6.12)式及式可化为因此，当 时，得由此得卫星的轨道是一个圆.即把地球半径，质量及引力常数的具体数值代入上式，并计算得即称为第一宇宙速度第一宇宙速度，此时当时，由(4.6.13)得即所求速度是第一宇宙速度的倍，即，称为第二宇宙速度第二宇宙速度，它的轨道是抛物线。是双曲线（一支）。当时，因，由(4.6.13)式可知这表明，当初速度小于第二宇宙速度时，卫星轨道是一个椭圆。当时，由(4.6.13)式可得因此，当初速度大于第二宇宙速度时，它的轨道由电路和机械装置组装成的永磁体扩音器模型如图4.7所示。4 4、扩音器振动模型、扩音器振动模型图4.7 一个时变电源电压 驱动着一个音圈能换器，能换器的转化系数为 ，通过能换器使扬声器的振动膜发生振动。由音圈组成的能换器，本质上是一个在永磁场内的自由运动。当变化的电流通过音圈时，音圈在永磁体的磁力和电流产生的磁力相互作用下进行运动。用 代表能换器与扬声器相互作用力，是能换器电阻，是能换器的感应系数.压定律，可得的微分方程组质量为的扬声器的振动是一个具有阻尼的弹簧系统振动。阻尼系数为，弹簧的弹性系数为，能换器转化系数与有关变量间的相互关系为这里是音圈两端的电压降；是音圈位移；代表音圈的速度，由牛顿第二定律及回路电(4.6.14)令，把(4.6.14)化为一阶微分方程组(4.6.15)方程组(4.6.15)的向量形式为，这里设 则矩阵 的特征方程为特征根为 ，其相应特征向量分别为因此，可求得齐次线形方程组的基解矩阵为为求(4.6.15)的通解，需求出(4.6.15)的一个特写成解，我们用待定系数法来求此特解。把因此，方程组(4.6.15)有特解将 代入方程 得比较上述方程两边和的系数得把(4.6.16)的消去得由此可求得和分别为的通解为方程这里 是任意常向量。进一步，可求得原方程中音圈的位移 和驱动能换器的电流 分别为5、人体含铅量问题、人体含铅量问题铅是如何进入体内的铅是如何进入体内的?消化道消化道 通常是铅吸收的主要途径，儿童由这一途径吸收的铅所占比例约为85。铅在肠道内通过主动转运和被动扩散两种方式由小肠吸收人血液。主动转运占吸收总量的80以上。主动转运依赖铅与肠粘膜上一种转运蛋白质结合，由转运蛋白作载体将铅转运人血液。被动扩散则是铅由肠腔向血液自然扩散。肠腔中铅浓度越高、血液中铅浓度越低，被动扩散的量就越大。呼吸道呼吸道 空气中含铅的灰尘经鼻孔吸人呼吸道后，一部分被鼻毛、气管纤毛和支气管纤毛、细支气管纤毛阻挡，最后以痰液的形式返回口腔，儿童将其咽人消化道，再由消化道吸收入血。另一部分特别微小的铅尘到达肺泡，沉积在肺泡里，然后由吞噬细胞等吸收，进入血液。皮肤皮肤 当我们接触有机铅的时候，铅会被皮肤直接吸收。铅中毒会引起神经、消化、循环系统紊乱，表现为贫血、腹痛、高血压等。血铅高组的儿童的总智商、操作智商、语言智商分别比低血铅组落后14分、14分和13分，而每升血液中的铅浓度上升100微克，儿童的身高将降低13厘米。低剂量的铅接触可以对人体的红细胞、肾脏、免疫系统、骨髓和中枢神经的功能产生不良影响，而所有这些影响发生前都可能没有明显的临床症状。铅中毒会影响婴幼儿最初站立、行走和说话的年龄，也可能引起孩子注意力涣散、记忆力减退、理解力降低及学习困难等。儿童体内的铅水平可分从血、骨、齿、尿、发检测到，儿童血铅水平分为5级。级：血铅值低于10微克分升，身体处于相对安全状态；级：血铅值为10微克19微克分升，属于轻度铅中毒。影响造血、神经传导和认知能力，儿童易出现头昏、烦躁、注意力涣散、多动、厌食、腹胀、轻度贫血；级：血铅值达到20微克44微克分升，为中度铅中毒。引起缺钙、缺锌、缺铁、免疫力低下，运动不协调，视力和听力受损，学习困难、智力下降，生长发育迟缓，贫血、腹绞痛、食癖、反应迟钝等；级：血铅值达到45微克69微克分升，就是重度铅中毒。可出现性格改变、易激怒、攻击性行为、运动失调、贫血、腹绞痛、高血压、心律失常和痴呆等；级：当血铅值大于70微克分升时，为极重度铅中毒，可导致脏器损害、铅性脑病、瘫痪、昏迷等。问题和模型问题和模型 铅是一种人体所必须的微量元素，但体内铅含量过多时就会引起铅中毒.铅是一种重金属元素，通过食物、饮料、空气进入人体，经过呼吸和消化系统后进入血液，再经过血液循环慢慢进入人体和骨头中.铅可以经过人体的排泄系统、通出出汗、剪头发、剪指甲排出体外.我们根据铅在人体内的变化情况将人体分为血液、组织、骨头3个仓室，铅在这三个仓室中的转化关系如图所示.x1(t)表示t时刻血液中的含铅量，x2(t)a02*x3(t).表示t时刻组织中的含铅量，x3(t)表示t时刻骨头中的含铅量.假设在单位时间内从环境经过消化、吸收系统进入血液的铅为L，从血液进入组织、骨头的铅分别为a31*x1(t)和a21*x1(t)，从组织、骨头再进入血液的铅分别为a13*x2(t)、a12*x2(t)，血液和骨头向外界排出的铅分别为a01*x1(t)和示意图示意图 (1)根据前面的假设建立人体血液、组织和骨头中含铅量所满足的微分方程组,并求其通解.(2)取时间单位为天，对一个志愿者在一种环境中的生活情况进行测定得a21=0.011，a12=0.012,a31=0.0039，a13=0.000035,a01=0.021,a02=0.016,L=49.3.假设该志愿者开始时体内的含铅量为0，求他体内血液、组织和骨头中含铅量随时间变化的关系，画出这些解曲线的图，并求当x+时这些函数的极限.(4)假设该志愿者在此环境中生活了365天后搬到了一个无铅的环境中去(不再有外界的铅进入体内，即L=0)，再讨论他体内血液、组织和骨头中含铅量随时间变化的关系，并画出0t1460时这些含铅量曲线的图形.简化图简化图模型模型模型 x=Ax+b,数据数据Lead Transfer Coefficients Lead Transfer Coefficients(Rabinowitz,et al.)Units:days-1a a2121=0.011=0.011a a1212=0.012=0.012from blood to tissue and backa a3131=0.0039=0.0039a a1313=0.000035=0.000035from blood to bone and backa a0101=0.021=0.021a a0202=0.016=0.016excretion from blood and tissue时间分为两段:0,365和365,1460#定义常数定义常数a01:=0.021;a02:=0.016;a12:=0.012;a13:=0.000035;a21:=0.011;a31:=0.0039;L:=49.3;#定义矩阵和向量定义矩阵和向量A:=matrix(3,3,-(a01+a21+a31),a12,a13,a21,-(a02+a12),0,a31,0,-a13);b:=matrix(3,1,L,0,0);#定义方程定义方程equn1:=diff(x1(t),t)=-(a01+a21+a31)*x1(t)+a12*x2(t)+a13*x3(t)+L;equn2:=diff(x2(t),t)=a21*x1(t)-(a02+a12)*x2(t);equn3:=diff(x3(t),t)=a31*x1(t)-a13*x3(t);#解初始值问题解初始值问题dsolve(equn1,equn2,equn3,x1(0)=0,x2(0)=0,x3(0)=0,x1(t),x2(t),x3(t);结果太复杂结果太复杂,利用迭加原理、特征值和特利用迭加原理、特征值和特征向量法征向量法#非齐次方程的特解非齐次方程的特解with(linalg):with(plots):xe:=evalm(-(inverse(A)&*b);#特征值和特征向量特征值和特征向量eigenvals(A);eigenvects(A);lambda:=-.3061796847e-4,-.1980315266e-1,-.4410122938e-1;v1:=.1124436946e-1,.442226656e-2,10.00746814;v2:=-.5975248926,-.8018660787,.1178839056;v3:=-.8256380196,.5640574390,.7307156328e-1;#特征向量组成矩阵特征向量组成矩阵 P:=augment(v1,v2,v3);#得到解的表达式得到解的表达式x1:=c1*P1,1*exp(lambda1*t)+c2*P1,2*exp(lambda2*t)+c3*P1,3*exp(lambda3*t)+xe1,1;x2:=c1*P2,1*exp(lambda1*t)+c2*P2,2*exp(lambda2*t)+c3*P2,3*exp(lambda3*t)+xe2,1;x3:=c1*P3,1*exp(lambda1*t)+c2*P3,2*exp(lambda2*t)+c3*P3,3*exp(lambda3*t)+xe3,1;#在解的表达式用在解的表达式用t=0t=0代入代入x10:=simplify(subs(t=0,x1);x20:=simplify(subs(t=0,x2);x30:=simplify(subs(t=0,x3);#求解常数求解常数solve(x10=0,x20=0,x30=0,c1,c2,c3);assign(%);#作图作图plot(x1,x2,x3,t=0.365,color=red,blue,green,thickness=2);plot1:=%:#重新显示解重新显示解x1;x2;x3;#求极限求极限limit(x1,t=infinity);limit(x2,t=infinity);limit(x3,t=infinity);#得到得到365365天后解的表达式天后解的表达式xx1:=cc1*P1,1*exp(lambda1*t)+cc2*P1,2*exp(lambda2*t)+cc3*P1,3*exp(lambda3*t);xx2:=cc1*P2,1*exp(lambda1*t)+cc