EN 560.601 - Spring 2019Homework #10Due on 12PM May 1st, 2019This problem set contains significant portion of coding. So please attachyour code and the result in your HW.1. Lotka-Volterra predator-prey modelConsider a continuous predator-prey model, where rabbit is the preyand fox is the predator. Let x(t) and y(t) be the population of therabbits and fox, respectively. The birth rate of the rabbits populationis proportional to its own population: αx(t). The death rate, on theother hand, is proportional both the rabbit and the fox populations:βx(t)y(t), meaning that the more the fox are, the more quickly therabbits will disappear. The reverse holds for the fox population: thebirth rates depends on the fox population as well as how many rabbitsare there as food: δx(t)y(t). The death rate of the fox population isproportional to its own population: γy(t). Thus the following LotkaVolterraequations are obtained:dxdt = αx βxydydt =γy + δxyIt is hard to find the analytic solutions for x(t) and y(t). However, aconserved “energy” E(x, y) can be defined for checking the consistencyof the numerical solutions:E(x, y) =δx + γ ln x βy + α ln yYou have done that in HW4.In this exercise, take α = γ = 0.5 and β = δ = 0.01.Use ode45 to solve the Lotka-Volterra equations and obtain x(t) andy(t) for t between t = 0 and t = 50, given at time t = 0, x(0) = 100and y(0) = 80. With the default settings, the number of time stepstaken should be 133.1ANSWERS: In the graph plot(t,x(t),t,y(t)), you should findthat the solutions x(t) and y(t) seems to be periodic and oscillatesomewhere between 10 and 120. These are the non-stationary steadystate solutions. By taking dx/dt = dy/dt = 0 in the Lotka-Volterraequations, one finds that the non-zero temporal means of x(t) andy(t) are γ/δ = α/β = 50. Therefore, the point (x, y) = (50, 50) isan “attractor” in the phase diagram plot(x(t),y(t)), and the tracex(t), y(t)should form a closed orbit around the attractor. However,if you make the phase diagram, you will see that the path is actuallynot closed due to numerical instabilities, i.e. the Lotka-Volterra equationsare stiff with respect to ode45. The numerical instabilities canbe revealed by checking the “constancy” of the energy function. If youplot E (x(t), y(t)) as a function of t, you will see the irregularity of theenergy evolution.* Define dE1 as the difference between the maximum and the minimumvalues of Ex(t), y(t)for the ode45 solutions.The stiffness of the system can be dealt with by using a better ODEsolver. Repeat the above calculations using the Trapezoidal RuleBackwardDifferentiation Formula (TR/BDF2) solver ode23tb. Thetrace in the phase diagram is “more” closed.* Define dE2 as the difference between the maximum and the minimumvalues of Ex(t), y(t)for the ode23tb solutions.Another way to deal with the numerical instability is to reduce thetolerance of the ODE solver. Use odeset to redefine the relative toleranceRelTol to 10?6 and calculate x(t), y(t)using ode45 with theoptions set by odeset.* Define dE3 as the difference between the maximum and the minimumvalues of Ex(t), y(t)for the ode45 solutions with the adjustedRelTol.Please plot everything as requested above and find dE1, dE2, dE3 .22. In class, we give an improved numerical method for solving ODE, calledHeun’s method. Heun’s method is a two step method. Here is thescheme,(a) Use the Taylor expansion to determine the local truncation error k. What is 代写EN 560.601作业、代做MATLAB程序设计作业、代写MATLAB课程作业、代做linear system作业the global truncation error Ek?(b) Consider the equation y0 = f(y) with f(y) = λy. Here λ is aconstant. Please derive the iteration relations for yi+1 and yi.Since we don’t want yi blows up to +∞, we can get the stabilitycondition for Heun’s method. Assume λ is a real number, pleasecompute the range of it where the scheme is stable as a functionof λ.(c) Implement this algorithm in MATLAB. Please test your algorithmwith the example presented in class,u000 + u2u0 + cos(t)u = 0, u(0) = 1, u0(0) = 1, u00(0) = 0in the time domain [0, 10] with dt = 0.02. Use ode45 to estimatethe true solution utrue in the same time domain. Find the L2error norm(u-utrue).3. Consider the boundary value problem (Note: this is S-L problem)with φn(4) = 0 and φn(4) = 0. Note φn is normalized to 1, i.e,2dx = 1. For this calculation, use x ∈ [4, 4] and choosexspan = 4 : 0.1 : 4. Find and plot the first three normalized eigenfunctionsand eigenvalues �n. Be sure to use a forward- and backwarddifferencingfor the boundary conditions.34. Consider the following linear system�This is a discrete version of Poisson’s equation: φ = ρ.Construct the A matrix above for a 99 × 99 system. Set k = 63 andconstruct a ρ (rho) vector according to the following rule:ρj = 2(1 cos(kπ/100)) sin(j kπ/100).(a) Please program the Jacobi and Gauss-Siedel methods for solvingthis problem. Start the iteration with an initial guess φ as a columnof ones. Continue to iterate until every term in the vector φconverges to within 10?4 of the same term in the previous iteration,i.e, norm(phi(:,i+1)-phi(:,i),inf)(b) In the last part, we find out the Jacobi method needs 5129 iterationsto reach the convergent condition. however, the Gauss-Seidelmethod only takes 2566 iterations to reach the same condition. Sothe Gauss-Seidel is almost twice faster than Jacobi method. Wewould like to show it by analyzing the eigenvalues of iterative matricesfor each method. Let’s remind the iterative systems for eachmethod. If the linear system of equations is Ax = b,Jacobi method: xk+1 = Db, where A = D T. (1)Gauss-Seidel method: xk+1 = Lb, where A = L U. (2)Define xt as the true solution of the linear system and the error asξk = xk xt. Show the iterations of xk will converge to the truesolution xtif the iterative systems are stable for both methods.Find the first 5 largest eigenvalues(in absolute value) of iterativematrices for each method in MATLAB. Are these two iterativesystems stable?45. This is the continuation of the last problem. There is another iterativemethod to solve the linear system of equations, called Successiveover-relaxation(SOR). This time the matrix A is decomposed into adiagonal matrix component D, and Strictly lower and upper triangularcomponents L and U:(This L matrix is different from the one in G-Smethod)A = D + L + U,�Its iteration relations isxk+1 = (D + ωL)(ωU + (ω 1)D)xk + (D + ωL)1ωbwhere ω is a constant, called relaxation factor. The range of ω is0 (a) Show when ω = 1, this method degenerates to Gauss-Seidel method.(b) Write down the iterative matrix for this method.(c) Please write a code to find the constant ω to minimize the spectralradius of this iterative matrix (the largest eigenvalue in absolutevalue) for the Poisson’s equation problem.(d) Use this optimal ω to implement SOR. How many iterations doyou need to reach the same convergent condition? How much thismethod faster than the Gauss-Seidel method?5转自：http://www.7daixie.com/2019050322005192.html
网友评论