代写EN 560.601、代做MATLAB程序设计、代写MATLAB、代做linear system

日期：2019-05-02 10:37

EN 560.601 - Spring 2019

Homework #10

Due on 12PM May 1st, 2019

This problem set contains significant portion of coding. So please attach

your code and the result in your HW.

1. Lotka-Volterra predator-prey model

Consider a continuous predator-prey model, where rabbit is the prey

and fox is the predator. Let x(t) and y(t) be the population of the

rabbits and fox, respectively. The birth rate of the rabbits population

is proportional to its own population: αx(t). The death rate, on the

other hand, is proportional both the rabbit and the fox populations:

βx(t)y(t), meaning that the more the fox are, the more quickly the

rabbits will disappear. The reverse holds for the fox population: the

birth rates depends on the fox population as well as how many rabbits

are there as food: δx(t)y(t). The death rate of the fox population is

proportional to its own population: γy(t). Thus the following LotkaVolterra

equations are obtained:

dx

dt = αx βxy

dy

dt =γy + δxy

It is hard to find the analytic solutions for x(t) and y(t). However, a

conserved “energy” E(x, y) can be defined for checking the consistency

of the numerical solutions:

E(x, y) =δx + γ ln x βy + α ln y

You 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) and

y(t) for t between t = 0 and t = 50, given at time t = 0, x(0) = 100

and y(0) = 80. With the default settings, the number of time steps

taken should be 133.

1

ANSWERS: In the graph plot(t,x(t),t,y(t)), you should find

that the solutions x(t) and y(t) seems to be periodic and oscillate

somewhere between 10 and 120. These are the non-stationary steady

state solutions. By taking dx/dt = dy/dt = 0 in the Lotka-Volterra

equations, one finds that the non-zero temporal means of x(t) and

y(t) are γ/δ = α/β = 50. Therefore, the point (x, y) = (50, 50) is

an “attractor” in the phase diagram plot(x(t),y(t)), and the trace

x(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 actually

not closed due to numerical instabilities, i.e. the Lotka-Volterra equations

are stiff with respect to ode45. The numerical instabilities can

be revealed by checking the “constancy” of the energy function. If you

plot E (x(t), y(t)) as a function of t, you will see the irregularity of the

energy evolution.

* Define dE1 as the difference between the maximum and the minimum

values of E

x(t), y(t)

for the ode45 solutions.

The stiffness of the system can be dealt with by using a better ODE

solver. Repeat the above calculations using the Trapezoidal RuleBackward

Differentiation Formula (TR/BDF2) solver ode23tb. The

trace in the phase diagram is “more” closed.

* Define dE2 as the difference between the maximum and the minimum

values of E

x(t), y(t)

for the ode23tb solutions.

Another way to deal with the numerical instability is to reduce the

tolerance of the ODE solver. Use odeset to redefine the relative tolerance

RelTol to 10?6 and calculate

x(t), y(t)

using ode45 with the

options set by odeset.

* Define dE3 as the difference between the maximum and the minimum

values of Ex(t), y(t)

for the ode45 solutions with the adjusted

RelTol.

Please plot everything as requested above and find dE1, dE2, dE3 .

2

2. In class, we give an improved numerical method for solving ODE, called

Heun’s method. Heun’s method is a two step method. Here is the

scheme,

(a) Use the Taylor expansion to determine the local truncation error k. What is the global truncation error Ek?

(b) Consider the equation y

0 = f(y) with f(y) = λy. Here λ is a

constant. Please derive the iteration relations for yi+1 and yi

.

Since we don’t want yi blows up to +∞, we can get the stability

condition for Heun’s method. Assume λ is a real number, please

compute the range of it where the scheme is stable as a function

of λ.

(c) Implement this algorithm in MATLAB. Please test your algorithm

with the example presented in class,

u

000 + u

2u

0 + cos(t)u = 0, u(0) = 1, u0

(0) = 1, u00(0) = 0

in the time domain [0, 10] with dt = 0.02. Use ode45 to estimate

the true solution utrue in the same time domain. Find the L2

error 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 choose

xspan = 4 : 0.1 : 4. Find and plot the first three normalized eigenfunctions

and eigenvalues n. Be sure to use a forward- and backwarddifferencing

for the boundary conditions.

3

4. 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 and

construct 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 solving

this problem. Start the iteration with an initial guess φ as a column

of 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)<=1e-4.

(b) In the last part, we find out the Jacobi method needs 5129 iterations

to reach the convergent condition. however, the Gauss-Seidel

method only takes 2566 iterations to reach the same condition. So

the Gauss-Seidel is almost twice faster than Jacobi method. We

would like to show it by analyzing the eigenvalues of iterative matrices

for each method. Let’s remind the iterative systems for each

method. If the linear system of equations is Ax = b,

Jacobi method: xk+1 = Db, where A = D T. (1)

Gauss-Seidel method: xk+1 = L

b, 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 true

solution xt

if the iterative systems are stable for both methods.

Find the first 5 largest eigenvalues(in absolute value) of iterative

matrices for each method in MATLAB. Are these two iterative

systems stable?

4

5. This is the continuation of the last problem. There is another iterative

method to solve the linear system of equations, called Successive

over-relaxation(SOR). This time the matrix A is decomposed into a

diagonal matrix component D, and Strictly lower and upper triangular

components L and U:(This L matrix is different from the one in G-S

method)

A = D + L + U,

Its iteration relations is

xk+1 = (D + ωL)

(ωU + (ω 1)D)xk + (D + ωL)

1ωb

where ω is a constant, called relaxation factor. The range of ω is

0 < ω < 2.

(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 spectral

radius of this iterative matrix (the largest eigenvalue in absolute

value) for the Poisson’s equation problem.

(d) Use this optimal ω to implement SOR. How many iterations do

you need to reach the same convergent condition? How much this

method faster than the Gauss-Seidel method?

5