2016暑期综合实训编程实验问题概要1讲解27231.pdf

《2016暑期综合实训编程实验问题概要1讲解27231.pdf》由会员分享，可在线阅读，更多相关《2016暑期综合实训编程实验问题概要1讲解27231.pdf（33页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、综合编程实验问题 第一章 常微分方程问题求解 第 1节 ODE45求解初值问题 1.基于数学建模中的传染病模型，应用 ODE45求解传染病 SIR问题，在同一图中画出 i(t)和 s(t)随 t 变化的曲线。对于不同的初始条件，在相平面中画出三条相轨线。2.基于数学建模中的种群竞争模型，种群依存模型和食饵捕食模型，应用ODE45求解模型，在时间空间和相平面上画出种群变化的图像，分析稳定点的稳定性。第 2节 编程计算 ODE 初值问题 1。编一个用 Euler方法解 atuTttutfu)(),(00 的程序，使之使用于任意右端函数f，任意步长h和任意区间,0Tt。用16/1,8/1,4/1h分

2、别计算初值问题 0666666.0151)0()4,0(,2utuuu 在结点)16,.,1,0(4/ii打印出问题的精确解（真解为)16/()(tteetu。计算近近似解、绝对误差、相对误差、先验误差界，分析输出结果（这与获得输出结果同样重要）2.编一个与上题同样要求的改进 Euler法的计算程序，1mu的初值用Euler方法提供，迭代步数s为输入参数。用它求解上题的问题，并将两个绍果加以比较。3。编一个程序用 Taylor级数法求解问题 1)0(10,uttuu 取 Taylor级数法的截断误差为)(21hO，即要用)(),.,(),()20(tututu的值 提示：可用一个简单的递推公式

3、来获得,.3,1),()(ntun 4。用四阶古典KuttaRunge方法（或其他精度不低于四阶的方法），对0 x时的标准正态分布函数 xdtexxt0,2121)(022 产生一张在5,0之间的80 个等距结点（即16/1h）处的函数值表。提示:寻找一个以)(x为解的初值问题。5。（一个“刚性”的微分方程）用四阶KuttaRunge方法解初值问题：0)0(30,1511102utttuu 取8/1h。每隔 8步打印出数值解与真解的值（tttu2/)(2），画出它们的大致图象，并对产生的结果作出解释。提示:当初值)0(u时，方程的真解为tteuttC220 6.分别用 Adams二步和四步外插

4、公式，用16/1h求解。1)0(30,17482utttuu 将计算结果与真解ttetut2/)(28进行比较，并对所产生的现象进行理论分析。7。用 Adams二步内插公式预测、Adams四步外插公式校正一次的预-校算法重新求解上题的方程、将结果与上题作比较并解释产生差异的原因。8。对（1.3）式所示的 Lotka-Volterra“弱肉强食模型，令 5,0,3,1,2,4tkdlekr，5,300yx，即 ttyxyxytyxyxtx0,5)0(,3)0(3)(24)(（l）取6/1h，用任过一种精度不低于三阶的办法求解，要求结果至少有三位有效数字。作出)(),(tytx的图像及y关于x的图

5、像。（2）对5.2,2,5.1,1)0(y解这同一个模型 分别画出y关于x的函数图象。（3）讨论所获得的结果并分析原因。提示:注意xy平面上的点（3，2）、它被称为平衡点)第 3节 常微分方程边值问题 1.调用函数 bvp4c 求解 MATLAB的的 5个例子，分析把高阶方程变为等价的一阶方程组的方法，剖析程序，总结编程求解过程。2.取64/1h和128/1h，计算以下两点边值问题的差分解，并与精确解比较（1）5.0)1(,1)0(10,)1(1)1(222uuxxuxdxud，精确解：xu11（2）3)(,2)0(0,sin322xxeuuxxeudxdudxud，精确解：u 精确解：xeu

6、xcos3。（3）)(,0)0(0,cos222uuxxxxudxduxdxud，精确解：xxusin2 并分析差分解与精确解的误差之所以会有些大有此小的原因。5 数值方法（英文版）习题和实验项目（ODE 数值解，9.1.3 习题）16考虑一阶微分方程)t(q)t(y)t(p)t(y 证明：一般解)t(y可用两个特殊积分求出。首先定义)t(F如下：dt)t(pe)t(F 然后，定义)t(y为 Cdt)t(q)t(F)t(F1)t(Y 提示：对乘积)t(y)t(F求导。17考虑放射物的衰减。如果)t(y是 t 时刻放射物的量，则)t(y将逐渐减少。实验表明，)t(y的变化率与未衰减物质的量成正比

7、。于是放射物衰减的初值问题为 0y)0(y,kyy (a)证明其解为kt0ey)t(y。(b)放射物质的半衰期是初始物质衰减一半所需的时间，14C的半衰期是 5730年。请给出求 t 时刻14C的量的公式)t(y。提示：求 k 使得0y5.0)5730(y.(c)分析一块木头后知，其中的14C的量是树木活着时的 0.712，该木头样本的年代有多久?(d)在某个时刻，一种放射物质的量为 10mg，23 s之后，该物质只剩1mgg。该物质的半衰期为多少 t 在习题 18 和习题 19 中，推导初值问题的方程并求解。18一个新的职业足球联赛的年度售票量计划以正比于 t 时刻的销售量和上限 3亿美元之

8、差的速度增长。假设最初的年售票量为 0美元，并且必须在 3年后达到 4000万美元(否则联赛取消)。基于这些假设，年销售量需要多久能达到 2200万美元?19一个新图书馆的内部容量为 5 百万立方英尺。通风系统以每分钟 4 5万立方英尺的速度引入新鲜空气。在通风系统打开之前，图书馆内部的二氧化碳和外面新鲜空气中的二氧化碳量分别为 0.4和 0.5。求通风系统打开 2小时之后图书馆中的二氧化碳百分比.9.2 欧拉方法 7 汪明当用欧拉方法求解b,a 上的初值问题 0y)a(y),t(fy0 时，结果为1M1kkh)t(f)b(y，它是逼近区间b,a 上)t(f的定积分的黎曼(Riemann)和。

9、8 说明欧拉方法不能求初值问题：0)0(y,y5.1)y,t(fy3/1 的近似解2/3t)t(y。证明你的结论，其中遇到了什么困难?9 能用欧拉方法求解0，3上的初值问题 0)0(y,y1y2 吗?提示：精确解为)ttan()t(y。p-7指数种群增长。某一种群以正比于当前数量的速度增长，且遵循O，5上的初值问题 5000)0(y,y02.0y (a)应用公式(10)，求出 y(5)的欧拉逼近，步长为 h=1，h=1/12和h=1/360.(b)(a)中当 h趋下 0时的极限是什么?p-8.一名跳伞运动员自飞机上跳下，降落伞打开之前的空气阻力正比于2/3v(v为速度)。设时间区间为O，6，向

10、下方向的微分方程为 0)0(v,v032.032v2/3 用欧拉方法和 h=0.05估计中 v(6)的值。p-9流行病模型。流行病的数学模型描述如下：设有 L 个成员的构成的群落，其中有 P个感染个体，Q为未感染个体。令)t(y)表示时刻 t 感染个体的数量。对于温和的疾病，如普通感冒，每个个体保持存活，流行病从感染者传播到未感染者。由于两组问有 PQ 种可能的 接触，)t(y的变化率正比于 PQ。故该问题可以描述为初值问题：0y)0(y),yL(kyy(a)用 L=25000，t=0.00003，h=0.2和初值条件，250)0(y，并用程序 9.1计算0，60上的欧拉近似解。(b)画出(a

11、)中的近似解。(c)通过求(a)中欧拉方法的纵坐标平均值来估计平均感染个体的数目。(d)通过用曲线拟合(a)中的数据，并用定理 1 10(积分均值定理)，估计平均感染个体的数目。P-10考虑一阶积分-常微分方程 t02d)(yy0001.0y25.0y3.1y (a)在区间O，20上，用欧拉方法和 h=0.2，y(0)=250以及梯形公式求方程的近似解。提示：欧拉方法的一般迭代公式为)d)(y.0001y0y25.0y3.1(hyykt0k2kkk1k 如果梯形公式用于逼近积分，则该表达式为)h(T.0001y0y25.0y3.1(hyykk2kkk1k 其中0)h(T0，且)yy(2h)h(

12、T)h(Tk1k1kk，99,.,1,0k (b)用初值y(O)=200和 y(O)=300重复(a)的计算。(c)在同一坐标系中画出(a)和(b)的近似解：13捕食者-被捕食者模型。非线性微分方程的一个例子是捕食者-被捕食者模型。设 x(t)和 y(t)分别表示兔子和狐狸在时刻 t 的数量，捕食者-被捕食者模型表明，)t(x和)t(y满足)t(Dy)t(y)t(Cx)t(y)t(y)t(Bx)t(Ax)t(x 一个典型的计算机模拟可使用系数 A=2，B=0.02，C=0.0002D=0.8 如果 (a)x(O)=3000只兔子，y(0)=120只狐狸 (b)x(0)=5000只兔子，y(0)

13、=100只狐狸 在区间0，5上用 M=50步和 h=0.2求解。6 case Study ODE Problems 1.The rate of change of the concentration of pollution in a lake is equal to the difference between the concentration of polluted water entering the lake and that leaving the lake.Assume that water containing a constant concentration of C kg/

14、km3 of pollutants enters the lake at a rate of 150 km3/year,and water leaves the lake at the same rate.Also assume that the volume of the lake remains constant at 5000km3.(a)Formulate a mathematical model to represent the rate of change of concentration of pollution in the lake.Find a mathematical s

15、olution.(b)If the initial concentration of pollution is 40 kg/km3,find the particular solution to the problem.(c)The fastest possible cleanup of the lake will occur if all pollution inflow ceases.This is represented by C=0.If all pollution into the lake was stopped immediately,how long would it take

16、 to reduce pollution to 50%of its current value?(d)Use the computer to graph your solution for the first 100 years after pollution stops.What happens to the concentration as time goes on?2.A projectile of mass 0.20kg is shot vertically upward with an initial velocity of 10 m/sec.It is then slowed do

17、wn due to the forces exerted by gravity and air resistance.(a)If the force due to air resistance equals 0.005 times the square of the projectiles instantaneous velocity acting in the opposite direction to the velocity,produce a mathematical model using an initial-value differential equation.Use velo

18、city as the dependent variable.Solve to find an expression for velocity in terms of time t.(b)Apply a fourth-order Runge-Kutta method with h=0.05 to estimate the projectiles instantaneous velocity for time t=0.05(0.05)1.00 sec.Validate your results using the exact solution from(a).3.An automobile sh

19、ock absorber coil spring system is designed to support 800lb,the portion of the automobiles weight it supports.The spring has a constant of 50 slugs/in.The effect of a bumpy road on the system can be described by the periodic function f(t)=300sin4t(in slug in/sec2),which acts upward on the tyre.The

20、system is initially in equilibrium at rest.(a)Assume that the automobiles shock absorber is so worn that it provides no effective damping force.Find a particular solution which describes the vertical displacement of the automobile over time.Use the computer to graph the particular solution for the f

21、irst ten seconds of motion.Describe the systems performance.(b)Now assume that the shock absorber is replaced.The new shock absorber exerts a damping force(in pounds)which is equal to 50 times the instantaneous vertical velocity of the system(in inches per second).Model this improved system with an

22、initial value problem.Solve it subject to the conditions described in part(a).Use the computer to graph the resulting equation for the first 10 seconds of the motion.Explain how the systems performance has improved.Is this system overdamped,underdamped or critically damped?4.A simple LRC electrical

23、circuit consists of a capacitor with a capacitance of 0.02 farads,a resistor with a resistance of 40 ohms and an inductor with an inductance of 8 henrys.The circuit is connected to a 24-volt battery.Initially there is no charge on the capacitor and no current in the circuit.Produce a mathematical mo

24、del which gives the charge on the capacitor for any time after the switch is closed.Find the charge on the capacitor after 1 sec and the current in the circuit after 1 sec.The LRC circuit is now connected to an alternating current source which applies a voltage E(t)=100cos 2t(in volts).There is no i

25、nitial charge on the capacitor or current in the circuit.(a)Find an equation which gives the charge on the capacitor for any time after the switch is closed.(b)Find the charge on the capacitor after 1 sec.(c)Find the current in the circuit after 1 sec.(d)Would a 5-ampere fuse have its capacity excee

26、ded in this circuit?(e)Use the computer to graph the transient response of the circuit(i.e.,the complementary function of the differential equation which models the circuit).(f)Use the computer to graph the general solution to the differential equation which models the circuit.Explain what happens t

27、o the transient response as time increases.Also,explain what happens to the steady state solution.5.Consider the following economic model.Let P be the price of a single item on the market and Q be the quantity of the item available on the market.Both P and Q are functions of time t.By considering pr

28、ice and quantity as two interacting species,the following mathematical model can be proposed:)QPc(QcdtdQPQcPcdtdP4321 where321c,c,c and4c are positive constants.Justify and discuss the adequacy of this model.(a)If 1c,10000c,1c321and 25c4,find the equilibrium points of this system.Classify each equil

29、ibrium point with respect to its stability.Give an explanation in cases where a point cannot be readily classified.(b)Use the computer to perform a graphical stability analysis to determine what will happen to the levels of P and Q as time increases.(c)Give an economic interpretation of the curves t

30、hat determine the equilibrium points.6.(a)For a simple RL circuit,Kirchhoffs voltage law requires that(if Ohms law holds)0RIdtdIL where L is inductance,R is resistance and I is current.Solve for I in the case L=R=2 and I(0)=0.01.Use both an analytical method and a numerical method.(b)In contrast to

31、part(a),real resistors may not always obey Ohms law.For example,the voltage drop may be nonlinear and the circuit dynamics may be described by a relationship such as 0RIII1dtdIL3refref where all other parameters are as defined in(a)and refI is a known reference current equal to 1.Solve for I as a fu

32、nction of time under the same conditions as specified in(a).7.Apart from inflow and outflow,another method by which mass can enter or leave a reactor is by a chemical reaction.For example,if the chemical decays,the reaction can sometimes be characterized as a first-order reaction,namely:Reaction=RVC

33、,where V=volume(m3),c=concentration(moles/m3)and R=reaction rate(min1),which can generally be interpreted as the fraction of the chemical which goes away per unit time.So,if R=0.1min-1,for example,then approximately 10%of the chemical in the reactor decays in one minute.On substituting the reaction

34、into the mass-balance equation,we have RVcFcFcdtdcVin where F=flow rate(m3/min).(a)Find the steady-state concentration of the reactor in the case where R=0.25 min1,cin=50mg/min3,F=10m3/min and V=200m3.(b)Repeat part(a),but compute the transient concentration response for the case 30m/mg20c.Validate

35、the results using Eulers numerical method fromt=0 to 30 min.8.Biomedical and environmental engineers must frequently predict the outcome of predator-prey or host-parasite relationships.A simple model of such interacts is provided by the following system of ODEs:212222211111xxgxddtdxxxdxgdtdx where 1

36、x and 2x are the numbers of hosts and parasites,respectively.The ds and gs are death and growth rates,respectively,where the subscript 1 refers to the host and the 2 to the parasite.Notice that the deaths of the host and the growth of the parasite are dependent on both x1 and x2.Use numerical method

37、s to compute values of)t(x1 and)t(x2 from t=0 to 10 for the following case:5.0d,01.0g,2.0d,1g,2)0(x,10)0(x221121 Use the computer to plot graphs of)t(x1 and)t(x2 against t.Interpret the results.9.A population of 1,000,000 people is subject to a disease which is seldom fatal and leaves the victim imm

38、une to future infections by this disease.Infection can only occur when a susceptible person comes into direct contact with an infectious person.The infectious period lasts approximately four weeks.Last week there were 25 new cases of the disease reported.This week there were 48 new cases.It is estim

39、ated that 25%of the population is immune due to previous exposure.(a)Develop a mathematical model as a discrete-time dynamical system.Hence find the eventual number of people who will become infected.(b)Estimate the maximum number of new cases in any one week.(c)Conduct a sensitivity analysis to inv

40、estigate the effect of any assumptions made in part(a)which were not supported by hard data.(d)Perform a sensitivity analysis for the number(25)of cases reported last week.It is thought by some that in early weeks the epidemic might be underreported.10.Consider a uniform beam of length l subject to

41、a linearly increasing distributed load lx0,l)x(W)x(W0 Assume the beam is hinged at the end x=0 and imbedded at the end x=l.By solving the governing ODE lx0,EI)x(Wdxyd44 where E is Youngs modulus of elasticity and I is the moment of inertia of the cross section about the neutral axis,show that the re

42、sulting deflection is given by)xlxl 2x(EIl120Wy43250 Taking the following parameter values:l=200in,E=29106lb/in2,I=725in4,W0=100lb/ft.use the computer to plot the elastic curve.Also use a numerical method to determine the point of maximum deflection,expressing your result in inches.第 2章 偏微分方程问题求解 1.

43、分析 MATLAB中的 PDEX1和 PDEX4的程序，总结 PDE 求解的一般步骤，并对 PDEX2至 PDEX5进行求解分析。2.设G是以原点为中心的单位正六边形的内部，用8/1h的正方形网格作剖分，用五点差分格式求方程：),(,0),(,1yxuGyxu 的数值解。3.对方程：),(,),(,0yxxuGyxu G的形状如图 3.11所示，其中曲线部分为单位圆的 1/4。取4/1h，求出所有内结点上的差分解。4.考虑图 3.12所示的的不规则区域上的热的分布问题：在两侧在底部在顶部在内部,0,100,00 xuuuu 决定相应的线性方程组并求出 10 个内部结点处的温度)10,.,2,1

44、(iui。5。对定解问题 0)0,(0),1(),0(,0(),1,0(,22xututuTtxxutu 若在)0,21(u处有一个扰动1021，取2/1,16/1rh，分别用古典显式格式和 Richardson格式计算 8层。（1）打印出第 8层 L各结点处的计算值；（2）预测继续算下去计算值的变化趋势；（3）分析上述趋桔产生的原因。6。用古典显格式求解定解问题：xkxututuTtxxutu)12sin()0,(0),1(),0(,0(),1,0(,22 分别取5/3,5/2r，取16/1h，计算 10-20层（1）对固定的k，比校5/2r和5/3r时的计算值的差值。（2）k分别取 1,2

45、,3,4,5，观测稳定和不稳定格式的计算值随初始函数变化的情况。7.改用 DuFort-Frankel格式（5.86）算出实习题 1在第 8 层的值，并与 Richardson格式的计算值作出比较。8.任选一种差分方法求自由振动问题),2(),0(sin)0,(,0)0,(0,2222tutuxxtuxutxxutu 的周期解，求出1t处一个周期的计算值。9.对定解问题 0),1(),0(1,2/1,12/1,0,)0,(,0(),1,0(,022tutuxxxxxuTtxxuaxutu 分别用表5-13的左偏显式和中心差显式。取16/1h，r分别为 3/5和 2/5,a分别为)3,2,1,0

46、(2nn计算 10 层，并分析所得到的计算结果，说说从山可获得什么规律性的东西。10.用16/1h，r分别为5/4和5/6的古典显格式计算 0),1(),0()12sin()0,(,0)0,(,0(),1,0(,2222tutuxkxtuxuTtxxutu)5,4,3,2,1(k，比较计算结果间的差别 5 Case Study PDE Problem 1.Write a computer program to determine the numerical solution of Laplaces equation 0yuxu2222 and Poissons equation fyuxu22

47、22 for a rectangular object of variable width and height.The object could have Dirichlet,Neumann or Cauchy boundary conditions.The value of f in Poissons equation should be assumed constant.Use this program to find the solution of the following problems:(a)A thin metal plate of dimension 2ft 2ft is

48、subjected to four heat sources which maintain the temperatures on its four edges as follows:u(x,0)=400oC,u(0,y)=200 oC,u(x,2)=50 oC,u(2,y)=100 oC.The flat sides of the plate are insulated so that no heat is transferred through these sides.Calculate the temperature profiles within the plate.(b)Perfec

49、t insulation is installed on two edges(right and top)of the plate of part(a).The other two edges are exposed to heat sources.This means that the set of Dirichlet and Neumann boundary conditions is 0|xu,0|yu,C200)y,0(u,C400)0,x(uy,22,xoo Calculate the temperature profiles within the plate and compare

50、 these with the results from part(a).(c)The thin metal plate of part(a)is made of an alloy which has a melting point ofC600o and a thermal conductivity of 15 Btu/(hour.ft.oC).The plate is subjected to an electric current which creates a uniform heat source within the plate.The amount of heat generat