綜合編程實驗問題第一章 常微分方程問題求解第1節(jié) ODE45求解初值問題1. 基于數(shù)學建模中的傳染病模型，應用ODE45求解傳染病SIR問題，在同一圖中畫出i(t)和s(t)隨t變化的曲線。對于不同的初始條件，在相平面中畫出三條相軌線。2. 基于數(shù)學建模中的種群競爭模型，種群依存模型和食餌捕食模型，應用ODE45求解模型，在時間空間和相平面上畫出種群變化的圖像，分析穩(wěn)定點的穩(wěn)定性。第2節(jié) 編程計算ODE初值問題1。編一個用Euler方法解的程序，使之使用于任意右端函數(shù)，任意步長和任意區(qū)間。用分別計算初值問題在結點打印出問題的精確解（真解為。計算近近似解、絕對誤差、相對誤差、先驗誤差界，分析輸出結果（這與獲得輸出結果同樣重要）2.編一個與上題同樣要求的改進Euler法的計算程序，的初值用Euler方法提供，迭代步數(shù)為輸入?yún)?shù)。用它求解上題的問題，并將兩個紹果加以比較。3。編一個程序用Taylor級數(shù)法求解問題取Taylor級數(shù)法的截斷誤差為，即要用的值 提示：可用一個簡單的遞推公式來獲得4。用四階古典方法（或其他精度不低于四階的方法），對時的標準正態(tài)分布函數(shù)產(chǎn)生一張在之間的80個等距結點（即）處的函數(shù)值表。提示:尋找一個以為解的初值問題。5。（一個“剛性”的微分方程）用四階方法解初值問題：取。每隔8步打印出數(shù)值解與真解的值（），畫出它們的大致圖象，并對產(chǎn)生的結果作出解釋。提示:當初值時，方程的真解為6.分別用Adams二步和四步外插公式，用求解。將計算結果與真解進行比較，并對所產(chǎn)生的現(xiàn)象進行理論分析。7。用Adams二步內(nèi)插公式預測、Adams四步外插公式校正一次的預-校算法重新求解上題的方程、將結果與上題作比較并解釋產(chǎn)生差異的原因。8。對（1.3）式所示的Lotka-Volterra“弱肉強食模型，令，即 （l）取，用任過一種精度不低于三階的辦法求解，要求結果至少有三位有效數(shù)字。作出的圖像及關于的圖像。（2）對解這同一個模型分別畫出關于的函數(shù)圖象。（3）討論所獲得的結果并分析原因。提示:注意平面上的點（3，2）、它被稱為平衡點)第3節(jié) 常微分方程邊值問題1.調(diào)用函數(shù)bvp4c 求解MATLAB的的5個例子，分析把高階方程變?yōu)榈葍r的一階方程組的方法，剖析程序，總結編程求解過程。2.取和，計算以下兩點邊值問題的差分解，并與精確解比較（1） ，精確解：（2） ，精確解：精確解：。（3） ，精確解：并分析差分解與精確解的誤差之所以會有些大有此小的原因。5 數(shù)值方法（英文版）習題和實驗項目（ODE 數(shù)值解，9.1.3 習題）16考慮一階微分方程 證明：一般解可用兩個特殊積分求出。首先定義如下：然后，定義為 提示：對乘積求導。17考慮放射物的衰減。如果是t時刻放射物的量，則將逐漸減少。實驗表明，的變化率與未衰減物質(zhì)的量成正比。于是放射物衰減的初值問題為 (a)證明其解為。 (b)放射物質(zhì)的半衰期是初始物質(zhì)衰減一半所需的時間，14C的半衰期是5730年。請給出求t時刻14C的量的公式。提示：求k使得. (c)分析一塊木頭后知，其中的14C的量是樹木活著時的0.712，該木頭樣本的年代有多久?(d)在某個時刻，一種放射物質(zhì)的量為10mg，23 s之后，該物質(zhì)只剩1mgg。該物質(zhì)的半衰期為多少t在習題18和習題19中，推導初值問題的方程并求解。18一個新的職業(yè)足球聯(lián)賽的年度售票量計劃以正比于t時刻的銷售量和上限3億美元之差的速度增長。假設最初的年售票量為0美元，并且必須在3年后達到4000萬美元(否則聯(lián)賽取消)。基于這些假設，年銷售量需要多久能達到2200萬美元?19一個新圖書館的內(nèi)部容量為5百萬立方英尺。通風系統(tǒng)以每分鐘45萬立方英尺的速度引入新鮮空氣。在通風系統(tǒng)打開之前，圖書館內(nèi)部的二氧化碳和外面新鮮空氣中的二氧化碳量分別為0. 4和0 .5。求通風系統(tǒng)打開2小時之后圖書館中的二氧化碳百分比.9.2 歐拉方法7汪明當用歐拉方法求解上的初值問題 時，結果為，它是逼近區(qū)間上的定積分的黎曼(Riemann)和。8說明歐拉方法不能求初值問題：的近似解。證明你的結論，其中遇到了什么困難?9能用歐拉方法求解0，3上的初值問題 嗎?提示：精確解為。p-7指數(shù)種群增長。某一種群以正比于當前數(shù)量的速度增長，且遵循O，5上的初值問題 (a)應用公式(10)，求出y(5)的歐拉逼近，步長為h=1，h=1/12和h=1/360. (b)(a)中當h趨下0時的極限是什么?p-8.一名跳傘運動員自飛機上跳下，降落傘打開之前的空氣阻力正比于(v為速度)。設時間區(qū)間為O，6，向下方向的微分方程為 用歐拉方法和h=0.05估計中v(6)的值。p-9流行病模型。流行病的數(shù)學模型描述如下：設有L個成員的構成的群落，其中有P個感染個體，Q為未感染個體。令)表示時刻t感染個體的數(shù)量。對于溫和的疾病，如普通感冒，每個個體保持存活，流行病從感染者傳播到未感染者。由于兩組問有PQ種可能的接觸，的變化率正比于PQ。故該問題可以描述為初值問題：(a)用L=25000，t=0.00003，h=0.2和初值條件，并用程序9.1計算0，60上的歐拉近似解。 (b)畫出(a)中的近似解。 (c)通過求(a)中歐拉方法的縱坐標平均值來估計平均感染個體的數(shù)目。 (d)通過用曲線擬合(a)中的數(shù)據(jù)，并用定理110(積分均值定理)，估計平均感染個體的數(shù)目。P-10考慮一階積分-常微分方程 (a)在區(qū)間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的數(shù)量，捕食者-被捕食者模型表明，和滿足一個典型的計算機模擬可使用系數(shù) A=2， B=0.02， C=0.0002D=0.8 如果 (a)x(O)=3000只兔子，y(0)=120只狐貍 (b)x(0)=5000只兔子，y(0)=100只狐貍 在區(qū)間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.設是以原點為中心的單位正六邊形的內(nèi)部，用的正方形網(wǎng)格作剖分，用五點差分格式求方程：的數(shù)值解。3.對方程：的形狀如圖3.11所示，其中曲線部分為單位圓的1/4。取，求出所有內(nèi)結點上的差分解。4.考慮圖3.12所示的的不規(guī)則區(qū)域上的熱的分布問題：決定相應的線性方程組并求出10個內(nèi)部結點處的溫度。5。對定解問題若在處有一個擾動，取，分別用古典顯式格式和Richardson格式計算8層。 （1）打印出第8層L各結點處的計算值； （2）預測繼續(xù)算下去計算值的變化趨勢； （3）分析上述趨桔產(chǎn)生的原因。6。用古典顯格式求解定解問題：分別取，取，計算10-20層（1）對固定的，比校和時的計算值的差值。（2）分別取1,2,3,4,5，觀測穩(wěn)定和不穩(wěn)定格式的計算值隨初始函數(shù)變化的情況。7.改用DuFort-Frankel格式（5.86）算出實習題1在第8層的值，并與Richardson格式的計算值作出比較。8.任選一種差分方法求自由振動問題的周期解，求出處一個周期的計算值。9.對定解問題分別用表5-13的左偏顯式和中心差顯式。取，分別為3/5和2/5,分別為計算10層，并分析所得到的計算結果，說說從山可獲得什么規(guī)律性的東西。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