讲解ODE留学生、辅导RKDP/IVP、讲解Matlab编程、辅导Matlab
- 首页 >> Matlab编程 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));
endLab 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 0sol = 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 0sol = 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 = tyLab 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.
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));
endLab 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
% 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
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 = tyLab 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.