代做ODE留学生、代写RKDP/IVP、代做Matlab编程、代写Matlab

日期：2019-03-05 09:06

Lab 4A: Solving ODE Initial-Value Problems

Due: Thursday, March 7th at 11:59PM

In this lab you will be numerically solving initial value ODEs. That is, differential equations of the form

y˙(t) = f

t, y(t)

on t0 ≤ t ≤ tfinal starting at t = t0 where we know y(t0) = y0 . You will use two numerical methods in

Matlab for solving ODEs.

ode45

The ode45 function in Matlab will compute a numerical solution of an IVP and works well for most

problems. The algorithm used is the Runge-Kutta Dormand-Prince (RKDP) method to produce fourth and

fifth order accurate solutions with only six function evaluations per step. This method works well for many

problems. For some problems, known as stiff problems, ode45 may not work as well. When ode45 fails or

does not produce an accurate solution the ode15s method should be tried.

ode15s

The ode15s method will produce better results when the problem is stiff. ode15s is a variable order solver

based on the numerical differentiation formulas (NDFs). Optionally, it uses the backward differentiation

formulas (BDFs, also known as Gear’s method) that are usually less efficient for non-stiff problems (taken

from the ode15s Matlab help page).

Stiff Problem

“A problem is stiff if ode15s mops the floor with ode45 on it." — L. F. Shampine

A stiff problem is one for which ode15s is more efficient than ode45. No, really. That is the practical way

to tell.

A stiff equation is a differential equation for which certain numerical methods for solving the equation are

numerically unstable, unless the step size is taken to be extremely small. It has proven difficult to formulate

a precise definition of stiffness, but the main idea is that the equation includes some terms that can lead to

rapid variation in the solution. (Wikipedia)

Rob does not like the Wikipedia definition. See [3].

Lab 4A: Solving ODE Initial-Value Problems 2

Converting Higher-Order ODE to a First-Order System

To convert a higher-order equation to a first-order system we use this one weird trick1

.

Example

y000 + 3y00 2 sin(y0) + y2 = 0 (1)

starting at y(0) = 0, y0

(0) = 1, and y

00(0) = 0. Because it is a third order differential equation, we introduce

a three-dimensional vector, say z, 3-dimensional because the equation was third order2

. Then put z1 = y,

z2 = y

0 and z3 = y

00 (up to one less than third order because we count like computer scientists: 0, 1, 2,

giving us three.) Then

z01 = z2 (because y0 = z2)z02 = z3 (because y00 = z3)z03 =3y

This trick is in Section 7.2 in Numerical Computing with Matlab. Also watch the Strang/Moler MIT

OCW videos!

Solving ODEs in Matlab

Example 1

Say we want to solve

on 0 ≤ t ≤ 1. Refer to [2] and [1] for an application of this ODE. This can be solved analytically, but most

ODEs cannot. The following Matlab code can be used to numerically solve this problem.

f = @(t, y) -sqrt(y);

sol = ode45(f, [0, 1], 1);

The sol variable in this case will be a special type of data structure used by the Matlab ODE solvers.

The following code can be used to interpolate the solution at 2019 points for t in [0, 1] .

t = linspace(0, 1, 2019);

y = deval(sol, t);

The solution can then be plotted.

plot(t, y);

Since the ODE was solved numerically we would like to know the error in the solution. Since we pretend

that we do not know the analytic solution to this ODE (and in general this would be right, because most of

the time we cannot find an analytic solution—if we could, why would we use numerics) we must look at the

backward error. The residual error (one kind of backward error) can be found from

1Which has already been covered in class, as a “companion matrix” for a polynomial.

2Three, count them, three.

Lab 4A: Solving ODE Initial-Value Problems 3

[y, dy] = deval(sol, t);

residual = dy - f(t, y)

which also computes the values of the derivative of the interpolant to the solution.

Now y is the exact solution to

y˙ = f(t, y) + residual

plot(t, residual, 'r.')

and we see this is small.

Example 2

Consider the system of equations

y˙1 = 2y1 + y2

y˙2 = y1 2y2 + y3

y˙3 = y2 2y3

with y(0) = [1, 0, 0]T on 0 ≤ t ≤ 5 . This system of ODEs can be solved in Matlab as follows.

% Set up the rhs of the system of ODEs

f = @(t, y) [-2*y(1) + y(2);

y(1) - 2*y(2) + y(3);

y(2) - 2*y(3)];

% Find the solution using ode15s

sol = ode15s(f, [0, 5], [1, 0, 0]');

% Interpolate the solution on 2019 points from 0 to 5

t = linspace(0, 5, 2019);

[y, dy] = deval(sol, t);

% Compute the residual

[~, m] = size(y);

residual = zeros(size(y));

for j=1:m

residual(:, j) = dy(:, j) - f(t(j), y(:,j));

end

An alternative way to solve this using matrices.

% Set up the rhs of the system of ODEs

A = [-2, 1, 0;

1, -2, 1;

0, 1, -2];

f = @(t, y) A*y;

% Find the solution using ode15s

sol = ode15s(f, [0, 5], [1, 0, 0]');

% Interpolate the solution on 2019 points from 0 to 5

t = linspace(0, 5, 2019);

[y, dy] = deval(sol, t);

% Compute the residual

[~, m] = size(y);

residual = zeros(size(y));

for j=1:m

residual(:, j) = dy(:, j) - f(t(j), y(:,j));

end

Lab 4A: Solving ODE Initial-Value Problems 4

Example 3

Consider the third-order equation

00(0) = 0. As was shown, this is equivalent to the system of first-order

equations

with z1(0) = 0, z2(0) = 1, and z3(0) = 0.

To solve this ODE in Matlab

% Set up the system of ODEs that is equivalent to the equation

f = @(t, z) [z(2);

z(3);

-3*z(3) + 2*sin(z(2)) - z(1).^2];

% Find the solution for 0<t<5

sol = ode45(f, [0, 5], [0, 1, 0]);

% Plot the solution, note that only the first row of z is plotted, this

% corresponds to z_1 which is the solution of the original third-order ODE.

% The second row of z is z_2 which is equal to y'.

t = linspace(0, 5, 2019);

z = deval(sol, t);

plot(t, z(1,:));

Example 4

Consider the first-order equation

y

0 = cos(πxy)

starting at 31 equally-spaced initial values on 0 ≤ x ≤ 5 and 0 ≤ t ≤ 5. To solve this ODE in Matlab

% Set up the system of ODEs that is equivalent to the equation

f = @(x, y) cos(pi*x*y);

% Find the solution for 0<t<5

sol = ode45(f, [0, 5], linspace(0, 5, 31));

xi = linspace(0, 5, 2019);

[z, dz] = deval(sol, xi);

% Plot the solution

plot(xi, z);

Lab 4A: Solving ODE Initial-Value Problems 5

Problem 1

In this problem you will be comparing the performance of ode45 and ode15s on the so-called ProtheroRobinson

test problem.

y

0 =λ(y cost) sin t , (2)

subject to y(0) = 100, for various values of λ.

Solve this system using both ode45 and ode15s for 0 ≤ t ≤ 10, for λ = 1, again for λ = 10, again for

λ = 100, and finally for λ = 1000 [Remember the minus sign in the equation.]. Make a semilogy plot of the

stepsizes from both methods for each of the four parameter values. [You may use a command something like

semilogy(sol.x(2:end), diff(sol.x), ’k.’ ), with whatever name you used, instead of ‘sol’.] Try

out the subplot function in Matlab to make side-by-side plots. It makes comparing the plots easier.

Similarly plot the residuals used by each method, again on a semilogy plot. Which method do you believe

did a better job at solving this problem? How do you know if this problem is stiff or not?

Problem 2

The Duffing equation is a non-linear second order differential equation that can be used to model some

damped and driven oscillators. See The Scholarpedia article on the Duffing Equation. Solve the Duffing

equation

x¨ + δx˙ + αx + βx3 = K cos(ωt)

with α = 1, β = 4, δ = 0.02, K = 8 and ω = 0.5 on 0 ≤ t ≤ 50 with initial condition x = ˙x = 1 .

You should make a plot of the solution and a separate plot of the residual.

In your report discuss the following as well as anything else you believe is important to solving this differential

equation.

Which ODE solver did you use for this problem and why? (feel free to try out some of the other ODE

solvers in Matlab such as ode23, ode113, etc.)

Is the problem stiff?

What does your residual plot mean?

You have not been given a theory of condition numbers for ODE (yet). Still, you have some tools in

your arsenal for estimating a condition number, using your judgement and by varying the problem

in some ways. Does this problem appear to be well-conditioned? Does the answer to that question

depend on the parameter values used here?

If your next meal depended on an accurate solution to this problem would you trust the solution you

found with Matlab? Why or why not?

Problem 3

This is an example of an ill-conditioned initial-value problem. Consider Airy’s differential equation

y

00 = ty

Lab 4A: Solving ODE Initial-Value Problems 6

subject to the initial condition

y(0) = Ai(0) = 13

2/3Γ(2/3)y0

(0) = Ai0

(0) =131/3Γ(1/3)

.

The reference solution is y(t) = Ai(t) which Matlab knows as airy. Matlab knows the Γ function by the

name gamma. The general solution of Airy’s equation is y = c1Ai(t) + c2Bi(t). For Bi(t) use airy(2, t) in

Matlab.

Solve Airy’s differential equation using ode45 on 0 ≤ t ≤ 2, and again on 0 ≤ t ≤ 10. Compare to airy(t).

Explain any discrepancy.

Submission Requirements

All files should be submitted on OWL on the Lab 4A submission. Only .m and .pdf files will be accepted

for grading. Submit all files needed to reproduce the results you found in the lab.

Your report should be written in a way such that someone with no programming experience can understand

what the problem was, how you solved it and what the results were.

References

[1] Robert M Corless and Julia Jankowski. Revisiting the discharge time of a cylindrical leaking bucket,

2016.

[2] Robert M Corless and Julia E Jankowski. Variations on a theme of Euler. SIAM Review, 58(4):775–792,

2016.

[3] Gustaf S derlind, Laurent Jay, and Manuel Calvo. Stiffness 1952–2012: Sixty years in search of a

definition. BIT Numerical Mathematics, 55(2):531–558, 2015.