1. No category

University of Washington
AMATH 301 Winter 2015
Instructor: Dr. King-Fai Li
Homework 4
Due: Friday Feb 27, 2015, 5 pm
Problem 1. First-order ordinary differential equations (ODEs)
Given the initial-value problem:
dy 2  2 yt
 2
,
dt
t 1
with
y  t  0  1 ,
we are to estimate the value of y  t  1 .
(a) Integrate the differential equation to t  1 using the forward Euler’s method with a time step
of h  0.1 . Save the resultant values  y  t  0  y  t  0.1  y  t  1  as a column
vector in A1.dat using the “save –ascii” command.
T
(b) Integrate the differential equation to t  1 using the 4th order Runge-Kutta method with a time
step of h  0.1 . Save the resultant values  y  t  0  y  t  0.1  y  t  1  as a column
vector in A2.dat using the “save –ascii” command.
T
(c) Find y  t  1 using the Matlab function ode45 with tspan=[0;1] and default settings
otherwise. Save the output of t and y  t  , in that order, as a 2-column matrix in A3.dat
using the “save –ascii” command. Note: If tspan has two elements only, then the
number of elements in the output t is automatically determined.
Problem 2. System of First-Order ODEs and Boundary-Value Problems
Let V represent the electrostatic potential between two concentric metal spheres of radii R1 and R2
(R1 < R2). The potential of the inner sphere is kept constant at V1 volts, and the potential of the
outer sphere is earthed to 0 volt. The potential in the region between the two spheres is governed
by Laplace’s equation, which, in this particular application, reduces to
d 2V 2 dV

 0,
dr 2 r dr
R1  r  R2 ,
V  R1   V1 ,
and
V  R2   0 .
(#)
Suppose R1 = 2 inches, R2 = 5 inches, and V1 = 110 volts. The purpose is to find V  r  for the inner
region between R1  r  R2 that satisfies the boundary conditions at R1 and R2. The system (#) is
called a boundary-value problem, which is different from the initial-value problems discussed in
dV  R1 
is specified instead of V  R2  . The ODE solvers in Matlab
the Video Lectures, where
dr
cannot handle boundary value problems. Below the shooting method will be used.
To start with, rephrase the above problem into an initial-value problem:
d 2V 2 dV

 0,
dr 2 r dr
R1  r  R2 ,
V  R1   V1 ,
and
dV
 R1    .
dr
(##)
where  is a constant to be determined such that the solution satisfies V  R2   0 . Define a vector
V 
 Y1  
Y     dV  . Then (##) can be transformed into a system of two first-order ODEs:

Y2  
 dr 
 Y 
d  Y1   2 
,
 2
dr Y2    Y2 
 r 
R1  r  R2 ,
Y1  R1   V1 ,
and
Y2  R1    .
(###)
In Matlab, (###) is equivalent to defining a function handle ode=@(r,Y)[ … ; … ];, where
dV
.
Y(1)  V and Y(2) 
dr
(a) Use ode45 to solve the system (###) with   100 volts/inch. The output Y in the
statement [r Y]=ode45(ode,…) has two columns: Y(:,1) corresponding to V  r  and
dV
 r  . Save the outputs of r and V  r  , in that order, as a 2dr
column matrix in A4.dat using the “save –ascii” command. Hint: V  r  R2  should
be negative.
Y(:,2) corresponding to
(b) Use ode45 to solve the system (###) with   90 volts/inch. Save the outputs of r and
V  r  , in that order, as a 2-column matrix in A5.dat using the “save –ascii” command.
Hint: V  r  R2  should be positive.
(c) The results obtained above suggest that the solution for  should lie between the interval
 100, 90 .
Use ode45 to calculate V  R2   Y(end,1) for   100, 99, , 90
 100 V  R2 

 100


 99 V  R2   99 
volts/inch. Save 
 as a 2-column matrix in A6.dat using the “save




 90 V R
 2   90 

–ascii” command.
(d) The final step is to estimate the solution for  using linear interpolation. Use interp1 to
T
finish this task. Treat V  R2   100 V  R2   99  V  R2   90  as x and [-100:-90]’


as y. Estimate the value of y when x=0 using y0=interp1(…). Save y0 in A7.dat using
the “save –ascii” command.
Problem 3. Stiff Systems
C. W. Gear proposed in 1971 a simple scalar example of stiffness:
dy
dG
   y  G  t   
,
dt
dt
  0,
dG
 1 , and
dt
where G  t  is a smooth, slowly varying external force. Let   41 , G  t   t ,
y  t  0   1 . We will integrate the system from t  0 to t  1 .
(a) Integrate the differential equation to t  1 using the forward Euler’s method with a time step
of h  0.05 . Save the resultant values  y  t  0  y  t  0.05  y  t  1  as a column
vector in A8.dat using the “save –ascii” command.
Lesson: The time series oscillates and amplifies because of the stiffness. A sufficiently small
h would make the amplification and oscillation disappear but then a lot of integration time
steps will be required.
T
(b) Integrate the differential equation to t  1 using the backward (implicit) Euler’s method with
a time step of h  0.05 . Save the resultant values  y  t  0  y  t  0.05  y  t  1  as
a column vector in A9.dat using the “save –ascii” command. Hint: Use fzero.
T
(c) Use ode23tb (implicit Trapezoidal Rule–Backward Difference) to solve the system with
tspan=[0;1] and default settings otherwise. Save the outputs of t and y  t  , in that order,
as a 2-column matrix in A10.dat using the “save –ascii” command. Hint: Number of
time steps < 50. Lesson: Just for fun, you can try ode45. The number of time steps taken
would be > 150 because ode45 should not be used for stiff systems.
Hint: Your solutions for (a) and (b) should look like
Problem 4. Lorenz model and Lyaponov constant
The Lorenz 1963 model for atmospheric convection is given by
dx
   y  x
dt
dy
 x  z  y
dt
dz
 xy   z
dt
  10
with
8
.
3
  28

In Matlab, define a function handle lorenz=@(t,Y)[ … ; … ; … ];, where Y(1)  x ,
Y(2)  y , and Y(3)  z .
(a) Use ode45 to solve the Lorenz model at t  0,0.01, 0.02,, 24.99, 25 with initial conditions
x  t  0   2 , y  t  0   3 , z  t  0   14 . Denote the corresponding outputs as t , xa  t  ,
ya  t  , and za  t  . Save t , xa  t  , ya  t  , and za  t  , in that order, as a 4-column matrix in
A11.dat using the “save –ascii” command. Hint: Define tspan=[0:0.01:25]. If
tspan has more than two elements, then the output t is the same as tspan. The appended
matrix [t Y] is what you need. Your solution should look like
plot3(Y(:,1),Y(:,2),Y(:,3)); view(2);
(b) To show the chaotic behavior of the system, one can integrate the Lorenz model again with
some tiny different initial conditions. Repeat (a) except z  t  0   14  10 9 . Denote the
corresponding outputs as t , xb  t  , yb  t  , and zb  t  . Save t , xb  t  , yb  t  , and zb  t  , in
that order, as a 4-column matrix in A12.dat using the “save –ascii” command. Note:
The output of t in (a) and (b) should be the same if the same tspan is set.
(c) The solutions  xa  t  , ya  t  , za  t   and  xb  t  , yb  t  , zb  t   will diverge away from each
other exponentially with time until a maximum distance is reached. Define the distance as
  t    xa  t   xb  t     ya  t   yb  t     za  t   zb  t   . Thus,   t   p2 exp  p1t  ,
where p1 is a Lyapunov exponent. Fit ln   t  with a straight line for 0  t  25 using
2
2
2
2
p=polyfit(…). The fitted slope is p1 . Save p as a row vector in A13.dat using the
“save –ascii” command. Note: p1 is the largest Lyapunov exponent of the Lorenz model,
which is close to 0.9. Your solution should look like
semilogy(t,Delta);
Note: After t = 25, ln   t  will level off because the distance cannot go beyond the size of the
stable orbit. In other words, one of the Lyapunov exponents must be zero.