2016暑期综合实训编程实验问题.doc

综合编程实验问题第一章 常微分方程问题求解第1节 ODE45求解初值问题1. 基于数学建模中的传染病模型，应用ODE45求解传染病SIR问题，在同一图中画出i(t)和s(t)随t变化的曲线。对于不同的初始条件，在相平面中画出三条相轨线。2. 基于数学建模中的种群竞争模型，种群依存模型和食饵捕食模型，应用ODE45求解模型，在时间空间和相平面上画出种群变化的图像，分析稳定点的稳定性。第2节 编程计算ODE初值问题1。编一个用Euler方法解的程序，使之使用于任意右端函数，任意步长和任意区间。用分别计算初值问题在结点打印出问题的精确解（真解为。计算近近似解、绝对误差、相对误差、先验误差界，分析输出结果（这与获得输出结果同样重要）2.编一个与上题同样要求的改进Euler法的计算程序，的初值用Euler方法提供，迭代步数为输入参数。用它求解上题的问题，并将两个绍果加以比较。3。编一个程序用Taylor级数法求解问题取Taylor级数法的截断误差为，即要用的值 提示：可用一个简单的递推公式来获得4。用四阶古典方法（或其他精度不低于四阶的方法），对时的标准正态分布函数产生一张在之间的80个等距结点（即）处的函数值表。提示:寻找一个以为解的初值问题。5。（一个“刚性”的微分方程）用四阶方法解初值问题：取。每隔8步打印出数值解与真解的值（），画出它们的大致图象，并对产生的结果作出解释。提示:当初值时，方程的真解为6.分别用Adams二步和四步外插公式，用求解。将计算结果与真解进行比较，并对所产生的现象进行理论分析。7。用Adams二步内插公式预测、Adams四步外插公式校正一次的预-校算法重新求解上题的方程、将结果与上题作比较并解释产生差异的原因。8。对（1.3）式所示的Lotka-Volterra“弱肉强食模型，令，即 （l）取，用任过一种精度不低于三阶的办法求解，要求结果至少有三位有效数字。作出的图像及关于的图像。（2）对解这同一个模型分别画出关于的函数图象。（3）讨论所获得的结果并分析原因。提示:注意平面上的点（3，2）、它被称为平衡点)第3节 常微分方程边值问题1.调用函数bvp4c 求解MATLAB的的5个例子，分析把高阶方程变为等价的一阶方程组的方法，剖析程序，总结编程求解过程。2.取和，计算以下两点边值问题的差分解，并与精确解比较（1） ，精确解：（2） ，精确解：精确解：。（3） ，精确解：并分析差分解与精确解的误差之所以会有些大有此小的原因。5 数值方法（英文版）习题和实验项目（ODE 数值解，9.1.3 习题）16考虑一阶微分方程 证明：一般解可用两个特殊积分求出。首先定义如下：然后，定义为 提示：对乘积求导。17考虑放射物的衰减。如果是t时刻放射物的量，则将逐渐减少。实验表明，的变化率与未衰减物质的量成正比。于是放射物衰减的初值问题为 (a)证明其解为。 (b)放射物质的半衰期是初始物质衰减一半所需的时间，14C的半衰期是5730年。请给出求t时刻14C的量的公式。提示：求k使得. (c)分析一块木头后知，其中的14C的量是树木活着时的0.712，该木头样本的年代有多久?(d)在某个时刻，一种放射物质的量为10mg，23 s之后，该物质只剩1mgg。该物质的半衰期为多少t在习题18和习题19中，推导初值问题的方程并求解。18一个新的职业足球联赛的年度售票量计划以正比于t时刻的销售量和上限3亿美元之差的速度增长。假设最初的年售票量为0美元，并且必须在3年后达到4000万美元(否则联赛取消)。基于这些假设，年销售量需要多久能达到2200万美元?19一个新图书馆的内部容量为5百万立方英尺。通风系统以每分钟45万立方英尺的速度引入新鲜空气。在通风系统打开之前，图书馆内部的二氧化碳和外面新鲜空气中的二氧化碳量分别为0. 4和0 .5。求通风系统打开2小时之后图书馆中的二氧化碳百分比.9.2 欧拉方法7汪明当用欧拉方法求解上的初值问题 时，结果为，它是逼近区间上的定积分的黎曼(Riemann)和。8说明欧拉方法不能求初值问题：的近似解。证明你的结论，其中遇到了什么困难?9能用欧拉方法求解0，3上的初值问题 吗?提示：精确解为。p-7指数种群增长。某一种群以正比于当前数量的速度增长，且遵循O，5上的初值问题 (a)应用公式(10)，求出y(5)的欧拉逼近，步长为h=1，h=1/12和h=1/360. (b)(a)中当h趋下0时的极限是什么?p-8.一名跳伞运动员自飞机上跳下，降落伞打开之前的空气阻力正比于(v为速度)。设时间区间为O，6，向下方向的微分方程为 用欧拉方法和h=0.05估计中v(6)的值。p-9流行病模型。流行病的数学模型描述如下：设有L个成员的构成的群落，其中有P个感染个体，Q为未感染个体。令)表示时刻t感染个体的数量。对于温和的疾病，如普通感冒，每个个体保持存活，流行病从感染者传播到未感染者。由于两组问有PQ种可能的接触，的变化率正比于PQ。故该问题可以描述为初值问题：(a)用L=25000，t=0.00003，h=0.2和初值条件，并用程序9.1计算0，60上的欧拉近似解。 (b)画出(a)中的近似解。 (c)通过求(a)中欧拉方法的纵坐标平均值来估计平均感染个体的数目。 (d)通过用曲线拟合(a)中的数据，并用定理110(积分均值定理)，估计平均感染个体的数目。P-10考虑一阶积分-常微分方程 (a)在区间O，20上，用欧拉方法和h=0. 2，y(0)=250以及梯形公式求方程的近似解。提示：欧拉方法的一般迭代公式为如果梯形公式用于逼近积分，则该表达式为其中，且， (b)用初值y(O)=200和y(O)=300重复(a)的计算。 (c)在同一坐标系中画出(a)和(b)的近似解：13捕食者-被捕食者模型。非线性微分方程的一个例子是捕食者-被捕食者模型。设x(t)和y(t)分别表示兔子和狐狸在时刻t的数量，捕食者-被捕食者模型表明，和满足一个典型的计算机模拟可使用系数 A=2， B=0.02， C=0.0002D=0.8 如果 (a)x(O)=3000只兔子，y(0)=120只狐狸 (b)x(0)=5000只兔子，y(0)=100只狐狸 在区间0，5上用M=50步和h=0.2求解。6 case Study ODE Problems1. 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 / km3 of pollutants enters the lake at a rate of 150 km3 / year , and water leavesthe 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 solution.(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 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 down 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 velocity 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 shock 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 functionf (t) =300sin4t (in slug in / sec2 ),which acts upward on the tyre. The 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 first 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 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 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 model 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 voltageE(t) =100cos 2t (in volts).There is no initial 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 exceeded 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 to 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 price and quantity as two interacting species, the following mathematical model can be proposed:where and are positive constants. Justify and discuss the adequacy of this model.(a) If and , find the equilibrium points of this system. Classify each equilibrium 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 that determine the equilibrium points.6. (a) For a simple RL circuit, Kirchhoffs voltage law requires that (if Ohms law holds)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 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 aswhere all other parameters are as defined in (a) and is a known reference current equal to 1. Solve for I as a function 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 ,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 into the mass-balance equation, we have 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 . Validate 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:where and 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 methods to compute values of and from t = 0 to 10 for the following case:Use the computer to plot graphs of and 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 immune 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 estimated that 25% of the population is immune due to previousexposure.(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 investigate 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 a linearly increasing distributed loadAssume the beam is hinged at the end x = 0 and imbedded at the end x = l .By solving the governing ODEwhere E is Youngs modulus of elasticity and I is the moment of inertia of the cross section about the neutral axis, show that the resulting deflection is given byTaking 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.分析MATLAB中的PDEX1和PDEX4的程序，总结PDE求解的一般步骤，并对PDEX2至PDEX5进行求解分析。2.设是以原点为中心的单位正六边形的内部，用的正方形网格作剖分，用五点差分格式求方程：的数值解。3.对方程：的形状如图3.11所示，其中曲线部分为单位圆的1/4。取，求出所有内结点上的差分解。4.考虑图3.12所示的的不规则区域上的热的分布问题：决定相应的线性方程组并求出10个内部结点处的温度。5。对定解问题若在处有一个扰动，取，分别用古典显式格式和Richardson格式计算8层。 （1）打印出第8层L各结点处的计算值； （2）预测继续算下去计算值的变化趋势； （3）分析上述趋桔产生的原因。6。用古典显格式求解定解问题：分别取，取，计算10-20层（1）对固定的，比校和时的计算值的差值。（2）分别取1,2,3,4,5，观测稳定和不稳定格式的计算值随初始函数变化的情况。7.改用DuFort-Frankel格式（5.86）算出实习题1在第8层的值，并与Richardson格式的计算值作出比较。8.任选一种差分方法求自由振动问题的周期解，求出处一个周期的计算值。9.对定解问题分别用表5-13的左偏显式和中心差显式。取，分别为3/5和2/5,分别为计算10层，并分析所得到的计算结果，说说从山可获得什么规律性的东西。10.用，分别为和的古典显格式计算，比较计算结果间的差别5 Case Study PDE Problem1. Write a computer program to determine the numerical solution of Laplaces equationand Poissons equationfor 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 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) Perfect 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 isCalculate the temperature profiles within the plate and compare 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 of 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 generated is Q = 100,000 Btu/(hour. ft3 ). The edges of the plate are in contact with heat sinks which maintain the temperature 50 oC on all four edges. Examine the temperature profiles within the plate to ascertain whether the alloy will begin to melt under these conditions.(d) Determine the optimum value of the overrelaxation parameter for the conditions used in part (a).2. Modify the computer program in Problem 1 to solve the three-dimensional problemApply this program to calculate the distribution of the dependent variable within a solid body which is subjected to the following boundary conditions:3. (a) Solve Laplaces equation with the following boundary conditionsDiscuss the results and determine the optimum value of the overrelaxation parameter for this problem.(b) Extend the computer program in (a) to include Robbins boundary conditions of the form:where u is the value of the dependent variable at the boundary and is a known value of the dependent variable in the fluid next to the boundary; k ,h and C are known constants.Apply this program to solve the following problem: The ambient temperature surrounding a house is 60 oF . The heating in the house has been turned off and so the internal temperature is also 60 oF at t = 0. The heating system is turned on and raises the internal temperature to 75 F at the rate of 5 F /hour.The ambient temperature remains at 60 F . The wall of the house is 0.4 ft thick and is made of material which has an average thermal diffusivity = 0.01ft2/hour and a thermal conductivity k = 0.2 Btu/(hour. ft2. F ). The heat transfer coefficient on the inside of the wall is h in =1.2Btu/(hour. ft2 . F) and the heat transfer coefficient on the outside is hout =2.0Btu/(hour. ft2 . F). Estimate how long it will take to reach a steadystatetemperature distribution across the wall.4. Develop the finite difference approximation of Ficks second law of diffusion in polar coordinates, namely,where c(r, ,z) represents the concentration and D the diffusivity. Hence write a computer program which can be used to solve the following problem:A wet cylinder of agar gel at 278 K with a uniform concentration of urea of 0.1 kg. mol/m3 has a diameter of 3cm and is 4cm long with flat parallel ends. The diffusivity is 4.51010m2/s . Calculate the concentration at the midpoint of the cylinder after 100 hours for the following cases if the cylinder is suddenly immersed in turbulent pure water:(a) Radial diffusion only.(b) Diffusion that occurs radially and axially.5. Consider a first-order chemical reaction being carried out under isothermal steady-state conditions, in a tubular-flow reactor. On the assumptions of laminar flow and negligible axial diffusion, the material balance equation iswhere= velocity of central stream line,=tube radius,reaction-velocity constant,radial diffusion constant,concentration of reactant,axial distance down the tube,radial distance from the centre.After defining the following dimensionless variables:the equation becomeswhere is the entering concentration of the reactant to the reactor.(a) Choose a set of appropriate boundary conditions for this problem and explain your choice.(b) What class of PDE is the above equation (hyperbolic, parabolic or elliptic)?(c) Set up the equation for numerical solution using finite differenceapproximations.(d) Does your choice of finite differences result in an explicit or implicit set of equations? Give details of the procedure for the solution of this set of equations.(e) Discuss stability considerations with respect to the method you have chosen.Figure 1. Stretched membrane fastened at the inside and outside boundaries. 6. A square membrane of side 12in (no bending or shear stresses), with a square hole of side 3in in the middle, is fastened at the inside and outside boundaries as shown in Figure 1. If a highly stretched membrane is subjected to a pressure p ,the PDE for the deflection u in the z -direction iswhere T is the tension (lb