Chapter 8 Ordinary Differential Equations:
We will now try to solve systems like this:
where f is sufficiently smooth with a variety of
numerical schemes.
8.1 Euler's Method:
This method is based on essentially tangent line
approximations. The iteration scheme generates
a sequence { yn } starting at y0 using
the following recursion equation:
| yn+1 |
= |
yn + h × f(tn,yn) |
| y0 |
= |
y0 |
where h is the step size we use for our
underlying partition of the time space giving
for appropriate indices.
8.1.1 The Matlab Implementation:
The basic code to implement the Euler method is quite
simple. Here it is:
FixedEuler.m
1 function [tvals,yvals] = FixedEuler(fname,y0,t0,tmax,h)
2 %
3 % fname fname is the name of the function f(t,y)
4 % h our choice of step size
5 % t0 starting point
6 % tmax end point
7 %
8 n = ceil((tmax-t0)/h)+1;
9 yvals = zeros(n,1);
10 tvals = zeros(n,1);
11 tvals(1) = t0;
12 for i=1:n-1
13 tvals(i+1) = tvals(i) + h;
14 end
15 if tvals(n) >= tmax
16 tvals(n) = tmax;
17 end
18 yvals(1) = y0;
19 for i=1:n-1
20 fval = feval(fname,tvals(i),yvals(i));
21 yvals(i+1) = yvals(i)+h*fval;
22 end
syntax highlighted by Code2HTML, v. 0.8.8b
8.1.2 The RunTime: Just Tables
We will test the Euler method code using two types of scripts;
one just does the Euler computation and the other will generate
a plot of the Euler approximation of the solution vs. the true
solution. Of course, this means you have to solve the ordinary
differential system for the solution by hand!
We will start by solving the system
for some choices of initial time t0 and initial state y0.
From basic differential equations, this system has
the solution
So we encode the right hand side function as func1.m and
the true solution as truefunc1.m with code
func1.m
1 function y = func1(t,x)
2 y = -x;
syntax highlighted by Code2HTML, v. 0.8.8b
and
truefunc1.m
1 function y = truefunc1(t0,y0,t)
2 y = y0*exp(t0-t);
syntax highlighted by Code2HTML, v. 0.8.8b
First, lets see what the Euler method generates as a table of values
via the script ShowEuler.m below
ShowEuler.m
1 while input('Another Step Size Choice? (1=yes, 0=no). ');
2 fname = input('Enter Function Name: ');
3 h = input('Input Step Size ');
4 t0 = input('Enter initial time ');
5 tmax = input('Enter final time: ');
6 y0 = input('Enter initial state: ');
7 freq = input('Input print frequency ');
8 s = ['Euler(' fname sprintf(',%6.3f,%6.3f,%6.3f,%6.3f)',y0,t0,tmax,h)];
9 disp([' tvals ' s])
10 disp(' ')
11 [tvals,yvals] = FixedEuler(fname,y0,t0,tmax,h);
12 N = size(tvals,1);
13 disp(sprintf(' %5.2f %20.16f',tvals(1),yvals(1)))
14 for m = freq+1:freq:N
15 disp(sprintf(' %5.2f %20.16f',tvals(m),yvals(m)))
16 end
17 disp(sprintf(' %5.2f %20.16f',tvals(N),yvals(N)))
18 end
syntax highlighted by Code2HTML, v. 0.8.8b
We then can test our Euler code by running this script.
Sample output is given below:
FixedEuler.prompts
1 >> ShowEuler
2 Another Step Size Choice? (1=yes, 0=no). 1
3 Enter Function Name: 'func1'
4 Input Step Size 0.4
5 Enter initial time 0.0
6 Enter final time: 5.0
7 Enter initial state: 1.0
8 Input print frequency 5
9 tvals Euler(func1, 1.000, 0.000, 5.000, 0.400)
10
11 0.00 1.0000000000000000
12 2.00 0.0777600000000000
13 4.00 0.0060466176000000
14 14
15 5.00 0.0013060694016000
16 Another Step Size Choice? (1=yes, 0=no). 0
syntax highlighted by Code2HTML, v. 0.8.8b
8.1.3 The RunTime: Plots!
Now let's generate some plots. We will use the following script
PlotEuler.m
1 while input('Another Step Size Choice? (1=yes, 0=no). ');
2 fname = input('Enter Function Name: ');
3 truename = input('Enter True Solution Name ');
4 h = input('Input Step Size ');
5 t0 = input('Enter Initial Time ');
6 tmax = input('Enter Final Time: ');
7 y0 = input('Enter Initial State: ');
8 s = ['Euler(' fname sprintf(',%6.3f,%6.3f,%6.3f,%6.3f)',y0,t0,tmax,h)];
9 disp([' tvals ' s])
10 disp(' ')
11 [tvals,yvals] = FixedEuler(fname,y0,t0,tmax,h);
12 N = size(tvals,1);
13 truevals = zeros(N,1);
14 for n=1:N
15 truevals(n) = feval(truename,t0,y0,tvals(n));
16 end
17 plot(tvals,truevals,tvals,yvals,'*');
18 print -deps eulerplot
19 end
syntax highlighted by Code2HTML, v. 0.8.8b
Here is the Matlab session:
PlotEuler.prompts
1 >> PlotEuler
2 Another Step Size Choice? (1=yes, 0=no). 1
3 Enter Function Name: 'func1'
4 Enter True Solution Name 'truefunc1'
5 Input Step Size 0.4
6 Enter Initial Time 1.0
7 Enter Final Time: 5.0
8 Enter Initial State: 1.0
9 tvals Euler(func1, 1.000, 1.000, 5.000, 0.400)
10
11 Another Step Size Choice? (1=yes, 0=no). 1
12 Enter Function Name: 'func1'
13 Enter True Solution Name 'truefunc1'
14 Input Step Size 0.1
15 Enter Initial Time 0.0
16 Enter Final Time: 5.0
17 Enter Initial State: 1.0
18 tvals Euler(func1, 1.000, 0.000, 5.000, 0.100)
19
20 Another Step Size Choice? (1=yes, 0=no). 0
syntax highlighted by Code2HTML, v. 0.8.8b
The plots for the case of stepsize h = 0.4 show a lot
of deviation from the true solution as we might expect since
we are just using a tangent line approximation.
The plot for the case of stepsize h = 0.1 is much better!
8.2 Runge-Kutta Methods:
Thes methods are based on multiple function
evaluations. The iteration scheme generates
a sequence { yn } starting at y0 using
the following recursion equation:
| yn+1 |
= |
yn + h × F(tn,yn,h,f) |
| y0 |
= |
y0 |
where h is the step size we use for our
underlying partition of the time space giving
for appropriate indices and F is some function of
the previous approximate solution, the step size and the
right hand side function f.
8.2.1 The Matlab Implementation:
The basic code to implement the Runge-Kutta methods is broken
into two pieces. The first one, RKstep.m implements
the evaluation of the next approximation solution at
point (tn,yn) given the old approximation at
(tn-1.yn-1).
Here is that code for Runge-Kutta codes of oders
one to four.
RKstep.m
1 function [tnew,ynew,fnew] = RKstep(fname,tc,yc,fc,h,k)
2 %
3 % fname the name of the right hand side function f(t,y)
4 % t is a scalar usually called time and
5 % y is a vector of size d
6 % yc approximate solution to y'(t) = f(t,y(t)) at t=tc
7 % fc f(tc,yc)
8 % h The time step
9 % k The order of the Runge-Kutta Method 1<= k <= 4
10 %
11 % tnew tc+h
12 % ynew approximate solution at tnew
13 % fnew f(tnew,ynew)
14 %
15 if k==1
16 k1 = h*fc;
17 ynew = yc+k1;
18 elseif k==2
19 k1 = h*fc;
20 k2 = h*feval(fname,tc+(h/2),yc+(k1/2));
21 ynew = yc + (k1+k2)/2;
22 elseif k==3
23 k1 = h*fc;
24 k2 = h*feval(fname,tc+(h/2),yc+(k1/2));
25 k3 = h*feval(fname,tc+h,yc-k1+2*k2);
26 ynew = yc+(k1+4*k2+k3)/6;
27 elseif k==4
28 k1 = h*fc;
29 k2 = h*feval(fname,tc+(h/2),yc+(k1/2));
30 k3 = h*feval(fname,tc+(h/2),yc+(k2/2));
31 k4 = h*feval(fname,tc+h,yc+k3);
32 ynew = yc+(k1+2*k2+2*k3+k4)/6;
33 else
34 disp(sprintf('The RK method %2d order is not allowed!',k));
35 end
36 tnew = tc+h;
37 fnew = feval(fname,tnew,ynew);
syntax highlighted by Code2HTML, v. 0.8.8b
Once the step is implemented, we solve the system
using the RK steps like this:
FixedRK.m
1 function [tvals,yvals] = FixedRK(fname,t0,y0,h,k,n)
2 %
3 % Gives approximate solution to
4 % y'(t) = f(t,y(t))
5 % y(t0) = y0
6 % using a kth order RK method
7 %
8 % t0 initial time
9 % y0 initial state
10 % h stepsize
11 % k RK order 1<= k <= 4
12 % n Number of steps to take
13 %
14 % tvals time values of form
15 % tvals(j) = t0 + (j-1)*h, 1 <= j <= n
16 % yvals approximate solution
17 % yvals(:j) = approximate solution at
18 % tvals(j), 1 <= j <= n
19 %
20 tc = t0;
21 yc = y0;
22 tvals = tc;
23 yvals = yc;
24 fc = feval(fname,tc,yc);
25 for j=1:n-1
26 [tc,yc,fc] = RKstep(fname,tc,yc,fc,h,k);
27 yvals = [yvals yc];
28 tvals = [tvals tc];
29 end
syntax highlighted by Code2HTML, v. 0.8.8b
8.2.2 The RunTime: Just Tables
We will test the RK methods code using two types of scripts;
one just does the RK computations and the other will generate
a plot of the RK approximation of the solution vs. the true
solution.
Again, we will solve the system
for some choices of initial time t0 and initial state y0.
From basic differential equations, this system has
the solution
We can then compare the RK solutions to the Euler solutions.
First, lets see what the RK method generates as a table of values
via the script ShowFixedRK.m below
ShowFixedRK.m
1 while input('Another Step Size Choice? (1=yes, 0=no). ');
2 fname = input('Enter Function Name: ');
3 h = input('Input Step Size ');
4 t0 = input('Enter initial time ');
5 tmax = input('Enter final time: ');
6 y0 = input('Enter initial state: ');
7 k = input('Enter RK Order in range [1,4]: ');
8 freq = input('Input print frequency ');
9 n = ceil((tmax-t0)/h)+1;
10 s = ['FixedRK(' fname sprintf(',%6.3f,%6.3f,%6.3f,%4d,%4d)',t0,y0,h,k,n)];
11 disp([' tvals ' s])
12 disp(' ')
13 [timevals,solutionvals] = FixedRK(fname,t0,y0,h,k,n);
14 tvals = timevals';
15 yvals = solutionvals';
16 disp(sprintf(' %5.2f %20.16f',tvals(1),yvals(1)))
17 for m = freq+1:freq:n
18 disp(sprintf(' %5.2f %20.16f',tvals(m),yvals(m)))
19 end
20 disp(sprintf(' %5.2f %20.16f',tvals(n),yvals(n)))
21 end
syntax highlighted by Code2HTML, v. 0.8.8b
We then can test our RK code by running this script.
Sample output is given below:
ShowFixedRK.prompts
1 >> ShowFixedRK
2 Another Step Size Choice? (1=yes, 0=no). 1
3 Enter Function Name: 'func1'
4 Input Step Size 0.4
5 Enter initial time 0.0
6 Enter final time: 5.0
7 Enter initial state: 1.0
8 Enter RK Order in range [1,4]: 2
9 Input print frequency 5
10 tvals FixedRK(func1, 0.000, 1.000, 0.400, 2, 14)
11
12 0.00 1.0000000000000000
13 2.00 0.1073741824000000
14 4.00 0.0115292150460685
15 5.20 0.0030223145490366
16 Another Step Size Choice? (1=yes, 0=no). 0
syntax highlighted by Code2HTML, v. 0.8.8b
8.2.3 The RunTime: Plots!
Now let's generate some plots. We will use the following script
PlotFixedRK.m
1 while input('Another Step Size Choice? (1=yes, 0=no). ');
2 fname = input('Enter Function Name: ');
3 truename = input('Enter True Solution Name ');
4 h = input('Input Step Size ');
5 t0 = input('Enter Initial Time ');
6 tmax = input('Enter Final Time: ');
7 y0 = input('Enter Initial State: ');
8 k = input('Enter RK Order in range [1,4]: ');
9 n = ceil((tmax-t0)/h)+1;
10 s = ['FixedRK(' fname sprintf(',%6.3f,%6.3f,%6.3f,%4d,%4d)',t0,y0,h,k,n)];
11 [timevals,solutionvals] = FixedRK(fname,t0,y0,h,k,n);
12 tvals = timevals';
13 yvals = solutionvals';
14 truevals = zeros(n,1);
15 for i=1:n
16 truevals(i) = feval(truename,t0,y0,tvals(i));
17 end
18 plot(tvals,truevals,tvals,yvals,'*');
19 print -deps fixedrkplot
20 end
syntax highlighted by Code2HTML, v. 0.8.8b
Here is the Matlab session:
PlotFixedRK.prompts
1 >> PlotFixedRK
2 Another Step Size Choice? (1=yes, 0=no). 11
3 Enter Function Name: 'func1'
4 Enter True Solution Name 'truefunc1'
5 Input Step Size 0.4
6 Enter Initial Time 0.0
7 Enter Final Time: 5.0
8 Enter Initial State: 1.0
9 Enter RK Order in range [1,4]: 2
10 Another Step Size Choice? (1=yes, 0=no). 1
11 Enter Function Name: 'func1'
12 Enter True Solution Name 'truefunc1'
13 Input Step Size 0.4
14 Enter Initial Time 0.0
15 Enter Final Time: 5.0
16 Enter Initial State: 1.0
17 Enter RK Order in range [1,4]: 3
18 Another Step Size Choice? (1=yes, 0=no). 1
19 Enter Function Name: 'func1'
20 Enter True Solution Name 'truefunc1'
21 Input Step Size 0.4
22 Enter Initial Time 0.0
23 Enter Final Time: 5.0
24 Enter Initial State: 1.0
25 Enter RK Order in range [1,4]: 4
26 Another Step Size Choice? (1=yes, 0=no). 0
syntax highlighted by Code2HTML, v. 0.8.8b
The plots for the case of stepsize h = 0.4 and RK order
2 are much better than the corresponding Euler plot.
The plots for the case of stepsize h = 0.4 and RK order
3 are even better:
The plots for the case of stepsize h = 0.4 and RK order
4 are the best yet!
8.3 The Adams Methods: Theory
We start with the fundamental system
Now integrate both sides and apply the initial condition to convert
this into an integral equation:
Now this means we can evaluate the above equation at the
partition points tn and tn+1:
Now subtract to obtain
This suggests that another way to approximate the solution of our
differential system is to replace the integrand f(t,y) by
an interpolating polynomial.
8.3.1 The Basics:
There are two types of
methods. We asume we are interested in calculating the new
point yn+1. There are two ways:
-
kth Order Adams - Bashford Explicit Methods
This method uses the previous k-1 values and the current
value yn as data points in the construction of our
interpolating polynomial. This means we use the domain data
| {(tn-k+1,yn-k+1), (tn-k+2,yn-k+2), ...,
(tn-1,yn-1, (tn,yn )} |
and function data
| {f(tn-k+1,yn-k+1), f(tn-k+2,yn-k+2), ...,
f(tn-1,yn-1), f(tn,yn)} |
to create an interpolant pk-1 with the pairs
| {(tn-k+1,fn-k+1), (tn-k+2,fn-k+2), ...,
(tn-1,fn-1, (tn,fn )} |
This method is called explicit because we will end up
being able to solve explicitly for yyn+1
with an equation of the form
g(yn+1 = g(fn, ..., fn-k+1).
- kth Order Adams - Moulton Implicit Methods
This method uses the previous k-2 values and the current
value yn as data points in the construction of our
interpolating polynomial. This means we use the domain data
| {(tn-k+2,yn-k+2), (tn-k+3,yn-k+3), ...,
(tn,yn, (tn+1,yn+1 )} |
and function data
| {f(tn-k+2,yn-k+2), f(tn-k+3,yn-k+3), ...,
f(tn,yn), f(tn+1,yn+1)} |
to create an interpolant pk-1 with the pairs
| {(tn-k+2,fn-k+2), (tn-k+3,fn-k+3), ...,
(tn,fn, (tn+1,fn+1 )} |
This method is called implicit because we end up with
an equation of the form g(yn+1,fn, ..., fn-k+2) = 0
which cannot be solved explicitly for yn+1 and instead
root finding tools must be used.
8.3.2 Adams-Bashford Methods:
We construct these explicit methods as follows:
| y(tn+1) |
= |
| y(tn) |: + |
ó
õ |
|
f(s,y(s)) ds |
|
| y(tn+1) |
approx
|
| y(tn) |: + |
ó
õ |
|
pk-1(s) ds |
|
where pk-1 is our kth order interpolant.
Consider table 8.1:
|
| Order |
Interpolant |
Data Points |
|
| 1 |
constant |
(tn,fn) |
| 2 |
linear |
(tn,fn), (tn-1,fn-1) |
| 3 |
quadratic |
(tn,fn), (tn-1,fn-1), (tn-2,fn-2) |
| 4 |
cubic |
(tn,fn), (tn-1,fn-1), (tn-2,fn-2), |
| |
|
(tn-3,fn-3) |
| 5 |
quartic |
(tn,fn), (tn-1,fn-1), (tn-2,fn-2), |
| |
|
(tn-3,yn-3), (tn-4,fn-4) |
Table 8.1: Adams - Bashford Methods
-
For the first order method, we approximate the integrand by
its value p0(t) = f(tn,yn) on the entire interval [tn,tn+1] to obtain
| yn+1 |
= |
yn + (tn+1 - tn) f(tn,yn) |
which is the usual Euler Method.
- For the second order method, we set
hn to be the difference tn+1 - tn and
Hence, we obtain,
| yn+1 - yn |
= |
| ó
õ |
|
|
æ
ç
ç
è |
fn-1 +
|
|
(s - tn-1)
|
ö
÷
÷
ø |
ds |
|
| |
= |
| hn fn-1 +
|
|
|
|
(tn+1 - tn-1)2
-
|
|
|
|
hn-12 |
|
| |
= |
| hn+1 fn-1 +
|
|
|
|
(hn + hn-1)2
-
|
|
|
|
hn-12 |
|
| |
··· |
|
| |
= |
|
|
( |
(hn + 2hn-1) fn - hn fn-1 |
) |
|
Note we can solve explicitly for the unknown yn+1.
Usually, we have uniform step sizes for our partition, so all of
the hn's are h. Using this we obtain
The derivation of the third, fourth and fifth (and higher!) methods
is quite similar but very messy. In Table 8.2,
we can see the results:
|
| Order |
Interpolant |
Equation |
|
| 1 |
constant |
yn+1 = yn +hfn |
| 2 |
linear |
yn+1 = yn +h/2(3fn - fn-1 |
| 3 |
quadratic |
yn+1 = yn + |
| |
|
h/12(23fn - 16 fn-1 + 5 fn-2) |
| 4 |
cubic |
yn+1 = yn + |
| |
|
h/24(55fn - 59fn-1 + 37fn-2 -9fn-3) |
| 5 |
quartic |
yn+1 = yn + |
| |
|
h/720(1901fn - 2774fn-1 + 2616fn-3 |
| |
|
- 1274fn-3 + 251fn-4) |
Table 8.2: Adams - Bashford Uniform Interpolants
8.3.3 Adams-Moulton Methods:
We construct these implicit methods as follows:
| y(tn+1) |
= |
| y(tn) |: + |
ó
õ |
|
f(s,y(s)) ds |
|
| y(tn+1) |
approx
|
| y(tn) |: + |
ó
õ |
|
pk-1(s) ds |
|
where pk-1 is our kth order interpolant.
Consider table 8.3:
|
| Order |
Interpolant |
Data Points |
|
| 1 |
constant |
(tn+1,fn+1) |
| 2 |
linear |
(tn+1,fn+1), (tn,fn) |
| 3 |
quadratic |
(tn+1,fn+1), (tn,fn), (tn-1,fn-1) |
| 4 |
cubic |
(tn+1,fn+1), (tn,fn), (tn-1,fn-1), |
| |
|
(tn-2,fn-2) |
| 5 |
quartic |
(tn+1,fn+1), (tn,fn), (tn-1,fn-1), |
| |
|
(tn-2,yn-2), (tn-3,fn-3) |
Table 8.3: Adams - Moulton Methods
-
For the first order method, we approximate the integrand by
its value p0(t) = f(tn+1,yn+1) on the entire interval
[tn,tn+1] to obtain
| yn+1 |
= |
yn + (tn+1 - tn) f(tn+1,yn+1) |
| |
= |
yn + hn f(tn+1,yn+1) |
Note that yn+1 appears on both sides so that we can't
solve for it. We instead must to a root finding on
| yn+1 - yn - hn f(tn+1,yn+1) |
= |
0 |
- For the second order method,
Hence, we obtain,
| yn+1 - yn |
= |
| ó
õ |
|
|
æ
ç
ç
è |
fn +
|
|
(s - tn)
|
ö
÷
÷
ø |
ds |
|
| |
= |
|
| |
= |
|
| |
= |
|
Note, if we write this out we see the implicit
nature of this calculation:
| yn+1 - yn
- |
|
|
( |
f(tn+1,yn+1) + f(tn,yn) |
) |
|
= |
0 |
For example, if f(t,y) = - t2y2, we would have
| yn+1 - yn
+ |
|
|
( |
tn+12 yn+12 + tn2 yn2 |
) |
|
= |
0 |
Usually, we have uniform step sizes for our partition, so all of
the hn's are h. Using this we obtain
The derivation of the third, fourth and fifth (and higher!) methods
is quite similar but very ugly. In Table 8.4,
we can see the results:
|
| Order |
Interpolant |
Equation |
|
| 1 |
constant |
yn+1 = yn +hfn+1 |
| 2 |
linear |
yn+1 = yn +h/2(fn + fn+1 |
| 3 |
quadratic |
yn+1 = yn + |
| |
|
h/12(5fn+1 + 8fn - fn-1) |
| 4 |
cubic |
yn+1 = yn + |
| |
|
h/24(9fn+1 + 19fn - 5fn-1 + 9fn-2) |
| 5 |
quartic |
yn+1 = yn + |
| |
|
h/720(251fn+1 + 646fn-1 - 264fn-2 |
| |
|
+ 106fn-3 - 19fn-4) |
Table 8.4: Adams - Moulton Uniform Interpolants
8.4 The Adams Bashford Methods in Matlab:
The Matlab code consist of the implementation and
several driver programs to illustrate their use
in solving a problem.
8.4.1 The Implementation:
The implementation of the Adams Bashford algorithm consists of
several parts:
-
ABstep.m: this code implements one AB step.
- StartAB.m: this code starts the kth order
AB method with the kth order RK method.
- FixedAB.m: this code computes the AB order k
solution to our given problem.
Here is the source code to the AB step implementation:
ABstep.m
1 function [tnew,ynew,fnew] = ABstep(fname,tc,yc,fvals,h,k)
2 %
3 % fname name of function f(t,y)
4 % t is scalar and y is column vector of size d
5 % yc approximate solution to y'(t) = f(t,y(t)) at tc
6 % fvals d by k matrix where fvals(:,i) is an approximation
7 % to f(t,y) at ti = tc +(1-i)h for 1 <= i <= k
8 % ie: at k-1 steps back, tc - (k-1)h to current time tc
9 % h time step
10 % k order of Adams Bashford method
11 %
12 % tnew tc+h
13 % ynew approximation at tnew
14 % fnew f(tnew,ynew)
15 %
16 if k ==1
17 ynew = yc+h*fvals;
18 elseif k==2
19 ynew = yc+(h/2)*(fvals*[3;-1]);
20 elseif k==3
21 ynew = yc+(h/12)*(fvals*[23;-16;5]);
22 elseif k==4
23 ynew = yc+(h/24)*(fvals*[55;-59;37;-9]);
24 elseif k==5
25 ynew = yc+(h/720)*(fvals*[1901;-2774;2616;-1274;251]);
26 else
27 disp(sprintf('The Adams Bashford method %2d order is not allowed!',k));
28 end
29 tnew = tc+h;
30 fnew = feval(fname,tnew,ynew);
syntax highlighted by Code2HTML, v. 0.8.8b
To start the AB methods, we use a RK method as follows:
StartAB.m
1 function [tvals,yvals,fvals] = StartAB(fname,t0,y0,h,k)
2 %
3 % Starts Adams Bashford Method of order k with corresponding
4 % order k Runge-Kutta method
5 %
6 % fname name of f(t,y)
7 % t0 initial time
8 % y0 initial state
9 % h step size
10 % k order of our methods
11 %
12 % tvals vector ti = t0 + (i-1)h, 1 <= i <= k
13 %
14 tc = t0;
15 yc = y0;
16 fc = feval(fname,tc,yc);
17 tvals = tc;
18 yvals = yc;
19 fvals = fc;
20
21 for j=1:k-1
22 [tc,yc,fc] = RKstep(fname,tc,yc,fc,h,k);
23 tvals = [tvals tc];
24 yvals = [yvals yc];
25 fvals = [fc fvals];
26 end
syntax highlighted by Code2HTML, v. 0.8.8b
Finally, we put it all together and compute the AB
order k solution.
FixedAB.m
1 function [tvals,yvals] = FixedAB(fname,t0,y0,h,k,n)
2 %
3 % fname name of f(t,y)
4 % t0 initial time
5 % y0 initial state
6 % h stepsize
7 % k order of our method
8 % n number of steps to take
9 %
10 % tvals vector of time values: ti = t0 + (i-1)h 1 <= i <= n
11 % yvals approximate solution at ti
12 %
13 [tvals,yvals,fvals] = StartAB(fname,t0,y0,h,k);
14 tc = tvals(k);
15 yc = yvals(:,k);
16 fc = fvals(:,k);
17
18 for j=k:n
19 %Take a step and then update
20 [tc,yc,fc] = ABstep(fname,tc,yc,fvals,h,k);
21 tvals = [tvals tc];
22 yvals = [yvals yc];
23 fvals = [fc fvals(:,1:k-1)];
24 end
syntax highlighted by Code2HTML, v. 0.8.8b
8.4.2 The RunTime: Just Tables
We can generate tables of AB output by running the
script ShowFixedAB.m. This script looks like:
ShowFixedAB.m
1 while input('Another Step Size Choice? (1=yes, 0=no). ');
2 fname = input('Enter Function Name: ');
3 h = input('Input Step Size ');
4 t0 = input('Enter initial time ');
5 tmax = input('Enter final time: ');
6 y0 = input('Enter initial state: ');
7 k = input('Enter RK Order in range [1,4]: ');
8 freq = input('Input print frequency ');
9 n = ceil((tmax-t0)/h)+1;
10 s = ['FixedAB(' fname sprintf(',%6.3f,%6.3f,%6.3f,%4d,%4d)',t0,y0,h,k,n)];
11 disp([' tvals ' s])
12 disp(' ')
13 [timevals,solutionvals] = FixedAB(fname,t0,y0,h,k,n);
14 tvals = timevals';
15 yvals = solutionvals';
16 disp(sprintf(' %5.2f %20.16f',tvals(1),yvals(1)))
17 for m = freq+1:freq:n
18 disp(sprintf(' %5.2f %20.16f',tvals(m),yvals(m)))
19 end
20 disp(sprintf(' %5.2f %20.16f',tvals(n),yvals(n)))
21 end
syntax highlighted by Code2HTML, v. 0.8.8b
and when we run the code, the Matlab session looks like this:
ShowFixedAB.prompts
1 >> ShowFixedAB
2 Another Step Size Choice? (1=yes, 0=no). 1
3 Enter Function Name: 'func1'
4 Input Step Size 0.4
5 Enter initial time 0.0
6 Enter final time: 5.0
7 Enter initial state: 1.0
8 Enter RK Order in range [1,4]: 2
9 Input print frequency 5
10 tvals FixedAB(func1, 0.000, 1.000, 0.400, 2, 14)
11
12 0.00 1.0000000000000000
13 2.00 0.1482240000000000
14 4.00 0.0231820697600000
15 5.20 0.0076120647270400
16 Another Step Size Choice? (1=yes, 0=no). 0
syntax highlighted by Code2HTML, v. 0.8.8b
8.4.3 The RunTime: Plots!
Now let's generate some plots. We will use the following script
PlotFixedAB.m
1 while input('Another Step Size Choice? (1=yes, 0=no). ');
2 fname = input('Enter Function Name: ');
3 truename = input('Enter True Solution Name ');
4 h = input('Input Step Size ');
5 t0 = input('Enter Initial Time ');
6 tmax = input('Enter Final Time: ');
7 y0 = input('Enter Initial State: ');
8 k = input('Enter AB Order in range [1,4]: ');
9 n = ceil((tmax-t0)/h)+1;
10 s = ['FixedAB(' fname sprintf(',%6.3f,%6.3f,%6.3f,%4d,%4d)',t0,y0,h,k,n)];
11 [timevals,solutionvals] = FixedAB(fname,t0,y0,h,k,n);
12 tvals = timevals';
13 yvals = solutionvals';
14 P = size(tvals,1);
15 truevals = zeros(P,1);
16 for i=1:P
17 truevals(i) = feval(truename,t0,y0,tvals(i));
18 end
19 plot(tvals,truevals,tvals,yvals,'*');
20 print -deps fixedabplot
21 end
syntax highlighted by Code2HTML, v. 0.8.8b
Here is the Matlab session:
PlotFixedAB.prompts
1 >> PlotFixedAB
2 Another Step Size Choice? (1=yes, 0=no). 1
3 Enter Function Name: 'func1'
4 Enter True Solution Name 'truefunc1'
5 Input Step Size 0.6
6 Enter Initial Time 0.0
7 Enter Final Time: 5.0
8 Enter Initial State: 1.0
9 Enter AB Order in range [1,4]: 2
10 Another Step Size Choice? (1=yes, 0=no). 1
11 Enter Function Name: 'func1'
12 Enter True Solution Name 'truefunc1'
13 Input Step Size 0.6
14 Enter Initial Time 0.0
15 Enter Final Time: 5.0
16 Enter Initial State: 1.0
17 Enter AB Order in range [1,4]: 4
18 Another Step Size Choice? (1=yes, 0=no). 1
19 Enter Function Name: 'func1'
20 Enter True Solution Name 'truefunc1'
21 Input Step Size 0.4
22 Enter Initial Time 0.0
23 Enter Final Time: 5.0
24 Enter Initial State: 1.0
25 Enter AB Order in range [1,4]: 3
26 Another Step Size Choice? (1=yes, 0=no). 0
syntax highlighted by Code2HTML, v. 0.8.8b
The plots for the case of stepsize h = 0.6 and AB order
2 are much better than the corresponding Euler plot.
The plots for the case of stepsize h = 0.6 and AB order
4 are even better, but there is still error:
We do much better by reducing the stepsize to 0.4
and order 3:
8.5 The Adams Moulton Methods in Matlab:
The Matlab code consist of the implementation and
several driver programs to illustrate their use
in solving a problem. We will work out the details
of what is called a predictor - corrector strategy here
and avoid the use of root finding tools to solve for
yk+1 in the AM schemes. Instead, we will
use the AB method to predict the value yk+1P
and then use that value on the right hand side of the
AM update and evaluate f(tk+1,yk+1P) in order
to compute the update yk+1 for the AM step. This
value is called the corrector value and is
labeled yk+1C.
8.5.1 The Implementation:
The implementation of the Adams Bashford algorithm consists of
several parts:
-
AMstep.m: this code implements one AB step.
- PCstep.m: this code finds the predictor
value from the AB step and the uses it in the
AM update equations to find the corrector value.
- FixedPC.m: this code computes the PC order k
solution to our given problem.
Here is the source code to the AM step implementation:
AMstep.m
1 function [tnew,ynew,fnew] = AMstep(fname,tc,yc,fvals,h,k)
2 %
3 % fname name of function f(t,y)
4 % t is scalar and y is column vector of size d
5 % yc approximate solution to y'(t) = f(t,y(t)) at tc
6 % fvals d by k matrix where fvals(:,i) is an approximation
7 % to f(t,y) at ti = tc +(2-i)h for 1 <= i <= k
8 % ie: at k-2 steps back, tc - (k-2)h to current time tc
9 % h time step
10 % k order of Adams Moulton method
11 %
12 % tnew tc+h
13 % ynew approximation at tnew
14 % fnew f(tnew,ynew)
15 %
16 if k ==1
17 ynew = yc+h*fvals;
18 elseif k==2
19 ynew = yc+(h/2)*(fvals*[1;1]);
20 elseif k==3
21 ynew = yc+(h/12)*(fvals*[5;8;-1]);
22 elseif k==4
23 ynew = yc+(h/24)*(fvals*[9;19;-5;1]);
24 elseif k==5
25 ynew = yc+(h/720)*(fvals*[251;646;-264;106;-19]);
26 else
27 disp(sprintf('The Adams Moulton method %2d order is not allowed!',k));
28 end
29 tnew = tc+h;
30 fnew = feval(fname,tnew,ynew);
syntax highlighted by Code2HTML, v. 0.8.8b
The Predictor - Corrector step is handled as follows:
PCstep.m
1 function [tnew,yPred,fPred,yCorr,fCorr] = ...
2 PCstep(fname,tc,yc,fvals,h,k)
3 %
4 % fname name of f(t,y)
5 % t is a scalar and y is a column vector of
6 % size d
7 % yc approximate solution to
8 % y'(t) = f(t,y(t)) at tc
9 % fvals d by k matrix where fvals(:i) is
10 % an approximation to f(t,y) at
11 % ti = to - (i-1)h for 1 <= i <= k
12 % h time step
13 % k Runge-Kutta order
14 %
15 % tnew tc + h
16 % yPred the predicted solution from AB at tnew
17 % fPred f(tnew,yPred)
18 % yCorr the corrected solution from AM at tnew
19 % fCorr f(tnew,yCorr)
20 %
21 [tnew,yPred,fPred] = ABstep(fname,tc,yc,fvals,h,k);
22 [tnew,yCorr,fCorr] = AMstep(fname,tc,yc,[fPred fvals(:,1:k-1)],h,k);
syntax highlighted by Code2HTML, v. 0.8.8b
Finally, we put it all together and compute the PC
order k solution.
FixedPC.m
1 function [tvals,yvals] = FixedPC(fname,t0,y0,h,k,n)
2 %
3 % fname name of f(t,y)
4 % t0 initial time
5 % y0 initial state
6 % h stepsize
7 % k order of our method
8 % n number of steps to take
9 %
10 % tvals vector of time values: ti = t0 + (i-1)h 1 <= i <= n
11 % yvals approximate solution at ti
12 %
13 [tvals,yvals,fvals] = StartAB(fname,t0,y0,h,k);
14 tc = tvals(k);
15 yc = yvals(:,k);
16 fc = fvals(:,k);
17
18 for j=k:n
19 %Take a step and then update
20 [tc,yPred,fPred,yc,fc] = PCstep(fname,tc,yc,fvals,h,k);
21 tvals = [tvals tc];
22 yvals = [yvals yc];
23 fvals = [fc fvals(:,1:k-1)];
24 end
syntax highlighted by Code2HTML, v. 0.8.8b
8.5.2 The RunTime: Just Tables
We can generate tables of PC output by running the
script ShowFixedPC.m. This script looks like:
ShowFixedPC.m
1 while input('Another Step Size Choice? (1=yes, 0=no). ');
2 fname = input('Enter Function Name: ');
3 h = input('Input Step Size ');
4 t0 = input('Enter initial time ');
5 tmax = input('Enter final time: ');
6 y0 = input('Enter initial state: ');
7 k = input('Enter PC Order in range [1,4]: ');
8 freq = input('Input print frequency ');
9 n = ceil((tmax-t0)/h)+1;
10 s = ['FixedPC(' fname sprintf(',%6.3f,%6.3f,%6.3f,%4d,%4d)',t0,y0,h,k,n)];
11 disp([' tvals ' s])
12 disp(' ')
13 [timevals,solutionvals] = FixedPC(fname,t0,y0,h,k,n);
14 tvals = timevals';
15 yvals = solutionvals';
16 disp(sprintf(' %5.2f %20.16f',tvals(1),yvals(1)))
17 for m = freq+1:freq:n
18 disp(sprintf(' %5.2f %20.16f',tvals(m),yvals(m)))
19 end
20 disp(sprintf(' %5.2f %20.16f',tvals(n),yvals(n)))
21 end
syntax highlighted by Code2HTML, v. 0.8.8b
and when we run the code, the Matlab session looks like this:
ShowFixedPC.prompts
1 >> ShowFixedPC
2 Another Step Size Choice? (1=yes, 0=no). 11
3 Enter Function Name: 'func1'
4 Input Step Size 0.5
5 Enter initial time 0.0
6 Enter final time: 5.0
7 Enter initial state: 1.0
8 Enter PC Order in range [1,4]: 2
9 Input print frequency 5
10 tvals FixedPC(func1, 0.000, 1.000, 0.500, 2, 11)
11
12 0.00 1.0000000000000000
13 2.50 0.0630731582641602
14 5.00 0.0041284898643426
15 5.00 0.0041284898643426
16 Another Step Size Choice? (1=yes, 0=no).
syntax highlighted by Code2HTML, v. 0.8.8b
8.5.3 The RunTime: Plots!
Now let's generate some plots. We will use the following script
PlotFixedPC.m
1 while input('Another Step Size Choice? (1=yes, 0=no). ');
2 fname = input('Enter Function Name: ');
3 truename = input('Enter True Solution Name ');
4 h = input('Input Step Size ');
5 t0 = input('Enter Initial Time ');
6 tmax = input('Enter Final Time: ');
7 y0 = input('Enter Initial State: ');
8 k = input('Enter PC Order in range [1,4]: ');
9 n = ceil((tmax-t0)/h)+1;
10 s = ['FixedPC(' fname sprintf(',%6.3f,%6.3f,%6.3f,%4d,%4d)',t0,y0,h,k,n)];
11 [timevals,solutionvals] = FixedPC(fname,t0,y0,h,k,n);
12 tvals = timevals';
13 yvals = solutionvals';
14 P = size(tvals,1);
15 truevals = zeros(P,1);
16 for i=1:P
17 truevals(i) = feval(truename,t0,y0,tvals(i));
18 end
19 plot(tvals,truevals,tvals,yvals,'*');
20 print -deps fixedpcplot
21 end
syntax highlighted by Code2HTML, v. 0.8.8b
Here is the Matlab session:
PlotFixedPC.prompts
1 >> PlotFixedPC
2 Another Step Size Choice? (1=yes, 0=no). 1
3 Enter Function Name: 'func1'
4 Enter True Solution Name 'truefunc1'
5 Input Step Size 0.6
6 Enter Initial Time 0.0
7 Enter Final Time: 5.0
8 Enter Initial State: 1.0
9 Enter PC Order in range [1,4]: 2
10 Another Step Size Choice? (1=yes, 0=no). 1
11 Enter Function Name: 'func1'
12 Enter True Solution Name 'truefunc1'
13 Input Step Size 0.6
14 Enter Initial Time 0.0
15 Enter Final Time: 5.0
16 Enter Initial State: 1.0
17 Enter PC Order in range [1,4]: 3
18 Another Step Size Choice? (1=yes, 0=no). 1
19 Enter Function Name: 'func1'
20 Enter True Solution Name 'truefunc1'
21 Input Step Size 1.0
22 Enter Initial Time 0.0
23 Enter Final Time: 5.0
24 Enter Initial State: 1.0
25 Enter PC Order in range [1,4]: 4
26 Another Step Size Choice? (1=yes, 0=no).
syntax highlighted by Code2HTML, v. 0.8.8b
The plots for the case of stepsize h = 0.6 and AB order
2 are just ok.
The plots for the case of stepsize h = 0.6 and AB order
3 are quite good: even though we use a huge step size!
If we go to a step size of 1.0 and order 4 we get:
8.6 Systems of Differential Equations:
Let's try to solve a system now. We will start with the
simple boundary value problem
on the interval [a,b]. This can be solved (although solving
a boundary value problem is not possible in general!) to give
| c |
= |
|
| d |
= |
y0 |
| f(t) |
= |
d + (t-a)(c+t2+at-2a2) |
We want to see how to solve this numerically.
8.7 How Do We Solve Systems?
We need to update our codes to handle systems.
8.7.1 Setting Up the Vector Functions:
The right hand side of our system is now a column vector
and the true solution is now the vector function
We will implement these vector functions in the Matlab
code vecfunc1.m and truevecfunc1.m respectively.
The right hand side is thus:
vecfunc1.m
1 function f = vecfunc1(t,y)
2 f = [y(2);6*t];
syntax highlighted by Code2HTML, v. 0.8.8b
while the true vector function is
truevecfunc1.m
1 function f = truefunc1(t0,t1,b,t)
2 c = ( b(2)-b(1) - t1^3 + t0^3 )/(t1-t0);
3 d = b(1)-t0^3;
4 f = [t^3+(t-t0)*(c-3*t0^2)+d;3*t^2+c-3*t0^2];
syntax highlighted by Code2HTML, v. 0.8.8b
8.7.2 Updating Our Solver Codes:
We will possibly need to update our codes
-
ShowFixedRK.m
- PlotFixedRK.m
Recall the ShowFixedRK.m script
ShowFixedRK.m
1 while input('Another Step Size Choice? (1=yes, 0=no). ');
2 fname = input('Enter Function Name: ');
3 h = input('Input Step Size ');
4 t0 = input('Enter initial time ');
5 tmax = input('Enter final time: ');
6 y0 = input('Enter initial state: ');
7 k = input('Enter RK Order in range [1,4]: ');
8 freq = input('Input print frequency ');
9 n = ceil((tmax-t0)/h)+1;
10 s = ['FixedRK(' fname sprintf(',%6.3f,%6.3f,%6.3f,%4d,%4d)',t0,y0,h,k,n)];
11 disp([' tvals ' s])
12 disp(' ')
13 [timevals,solutionvals] = FixedRK(fname,t0,y0,h,k,n);
14 tvals = timevals';
15 yvals = solutionvals';
16 disp(sprintf(' %5.2f %20.16f',tvals(1),yvals(1)))
17 for m = freq+1:freq:n
18 disp(sprintf(' %5.2f %20.16f',tvals(m),yvals(m)))
19 end
20 disp(sprintf(' %5.2f %20.16f',tvals(n),yvals(n)))
21 end
syntax highlighted by Code2HTML, v. 0.8.8b
We now run this script using the intial data
This generates the Matlab session
ShowFixedVecRK.prompts
1 >> ShowFixedRK
2 Another Step Size Choice? (1=yes, 0=no). 1
3 Enter Function Name: 'vecfunc1'
4 Input Step Size 0.2
5 Enter initial time 0.0
6 Enter final time: 1.0
7 Enter initial state: [1;2]
8 Enter RK Order in range [1,4]: 4
9 Input print frequency 5
10 tvals FixedRK(vecfunc1, 0.000, 1.000, 2.000,2.000000e-01, 4), 6.000,
11
12 0.00 1.0000000000000000
13 1.00 4.0000000000000000
14 1.00 4.0000000000000000
15 Another Step Size Choice? (1=yes, 0=no).
syntax highlighted by Code2HTML, v. 0.8.8b
and everything looks fine except that our header
s = ['FixedRK(' fname sprintf(',%6.3f,%6.3f,%6.3f,%4d,%4d)',t0,y0,h,k,n)];
will no longer work as we are trying to print out the vector
y0 into a single float data field! This could be fixed, but
note the fix would require that we know the size of the vector to
set up the print statement. Similar problems arise
when we try to use the plotting script from before.
Recall our original plotting script:
PlotFixedRK.m
1 while input('Another Step Size Choice? (1=yes, 0=no). ');
2 fname = input('Enter Function Name: ');
3 truename = input('Enter True Solution Name ');
4 h = input('Input Step Size ');
5 t0 = input('Enter Initial Time ');
6 tmax = input('Enter Final Time: ');
7 y0 = input('Enter Initial State: ');
8 k = input('Enter RK Order in range [1,4]: ');
9 n = ceil((tmax-t0)/h)+1;
10 s = ['FixedRK(' fname sprintf(',%6.3f,%6.3f,%6.3f,%4d,%4d)',t0,y0,h,k,n)];
11 [timevals,solutionvals] = FixedRK(fname,t0,y0,h,k,n);
12 tvals = timevals';
13 yvals = solutionvals';
14 truevals = zeros(n,1);
15 for i=1:n
16 truevals(i) = feval(truename,t0,y0,tvals(i));
17 end
18 plot(tvals,truevals,tvals,yvals,'*');
19 print -deps fixedrkplot
20 end
syntax highlighted by Code2HTML, v. 0.8.8b
To get this to work, we need to do several things:
-
We need to decide which component of our vectors
to plot.
- We need to worry more about row vs. column orientations
of our vectors in our calculations.
- We need to distinguish between the data vector
which controls the solution for the boundary value
problem and the intitial data vector which controls
the numerical solution of the ODE.
This is our altered code:
PlotVecFixedRK.m
1 while input('Do you want to Solve This Problem? (1=yes, 0=no). ');
2 fname = input('Enter Function Name: ');
3 truename = input('Enter True Solution Name: ');
4 dim = input('Enter Dimension of State: ');
5 h = input('Input Step Size: ');
6 t0 = input('Enter Initial Time: ');
7 tmax = input('Enter Final Time: ');
8 y0 = input('Enter Initial Data: ');
9 b = input('Enter Boundary Data: [value at to;value at t1]');
10 index = input('Enter Which Component to Plot: ');
11 k = input('Enter RK Order in range [1,4]: ');
12 n = ceil((tmax-t0)/h)+1;
13 s = ['FixedRK(' fname sprintf(',%6.3f,%6.3f,%6.3f,%4d,%4d)',t0,y0,h,k,n)];
14 [timevals,solutionvals] = FixedRK(fname,t0,y0,h,k,n);
15 tvals = timevals';
16 yvals = solutionvals';
17 truevals = zeros(n,dim);
18 for i=1:n
19 w = feval(truename,t0,tmax,b,tvals(i))';
20 truevals(i,:) = w;
21 end
22 plot(tvals,truevals(:,index),tvals,yvals(:,index),'*');
23 print -deps fixedvecrkplot
24 end
syntax highlighted by Code2HTML, v. 0.8.8b
Here is a run-time test: we will try to zero in on the
right initial condition to use to solve the boundary
value problem:
PlotVecFixedRK.prompts
1 >> PlotVecFixedRK
2 Do you want to Solve This Problem? (1=yes, 0=no). 1
3 Enter Function Name: 'vecfunc1'
4 Enter True Solution Name: 'truevecfunc1'
5 Enter Dimension of State: 2
6 Input Step Size: 0.2
7 Enter Initial Time: 0.0
8 Enter Final Time: 1.0
9 Enter Initial Data: [1;2]
10 Enter Boundary Data: [value at to;value at t1][1;2]
11 Enter Which Component to Plot: 1
12 Enter RK Order in range [1,4]: 4
13 Do you want to Solve This Problem? (1=yes, 0=no). 1
14 Enter Function Name: 'vecfunc1'
15 Enter True Solution Name: 'truevecfunc1'
16 Enter Dimension of State: 2
17 Input Step Size: 0.2
18 Enter Initial Time: 0.0
19 Enter Final Time: 1.0
20 Enter Initial Data: [1;0.05]
21 Enter Boundary Data: [value at to;value at t1][1;2]
22 Enter Which Component to Plot: 1
23 Enter RK Order in range [1,4]: 4
syntax highlighted by Code2HTML, v. 0.8.8b
with associated plots
We use initial data
Next, we use initial data
Clearyly, initial data
does the job.
8.8 Solving The BVP:
We solve the BVP by using what is called the shooting method.
Here is some Matlab code:
PlotShooting.m
1 while input('Do you want to Solve This Problem? (1=yes, 0=no). ');
2 fname = input('Enter Function Name: ');
3 truename = input('Enter True Solution Name: ');
4 dim = input('Enter Dimension of State: ');
5 h = input('Input Step Size: ');
6 t0 = input('Enter Initial Time: ');
7 tmax = input('Enter Final Time: ');
8 y0 = input('Enter Initial Data: ');
9 b = input('Enter Boundary Data: [value at to;value at t1]');
10 index = input('Enter Which Component to Plot: ');
11 k = input('Enter RK Order in range [1,4]: ');
12 n = ceil((tmax-t0)/h)+1;
13 s = ['FixedRK(' fname sprintf(',%6.3f,%6.3f,%6.3f,%4d,%4d)',t0,y0,h,k,n)];
14 [timevals,solutionvals] = FixedRK(fname,t0,y0,h,k,n);
15 tvals = timevals';
16 yvals = solutionvals';
17 error1 = (b(2) - yvals(n,1))^2;
18 y0new = input('Enter New Initial Data: ');
19 [timevals,solutionvals] = FixedRK(fname,t0,y0new,h,k,n);
20 tvals = timevals';
21 yvalsnew = solutionvals';
22 error2 = (b(2) - yvalsnew(n,1))^2;
23 c = y0(1)-error1*(y0new(2)-y0(2))/(error2-error1);
24 y0update = [b(1);c];
25 [timevals,solutionvals] = FixedRK(fname,t0,y0update,h,k,n);
26 tvals = timevals';
27 yvalsupdate = solutionvals';
28 error3 = (b(2) - yvalsupdate(n,1))^2;
29 truevals = zeros(n,dim);
30 for i=1:n
31 w = feval(truename,t0,tmax,b,tvals(i))';
32 truevals(i,:) = w;
33 end
34 disp(sprintf('error1 = %12.6f error2 = %12.6f error 3 = %12.6f',...
35 error1,error2,error3));
36 disp(sprintf('updated initial condition is c = %12.6f',c));
37 plot(tvals,truevals(:,index),tvals,yvalsupdate(:,index),'*');
38 print -deps shootingplot
39 end
syntax highlighted by Code2HTML, v. 0.8.8b
The run-time session is
PlotShooting.prompts
1 >> PlotShooting
2 Do you want to Solve This Problem? (1=yes, 0=no). 1
3 Enter Function Name: 'vecfunc1'
4 Enter True Solution Name: 'truevecfunc1'
5 Enter Dimension of State: 2
6 Input Step Size: 0.2
7 Enter Initial Time: 0.0
8 Enter Final Time: 1.0
9 Enter Initial Data: [1;2]
10 Enter Boundary Data: [value at to;value at t1][1;2]
11 Enter Which Component to Plot: 1
12 Enter RK Order in range [1,4]: 4
13 Enter New Initial Data: [1;1]
14 error1 = 4.000000 error2 = 1.000000 error 3 = 0.111111
15 updated initial condition is c = -0.333333
16 Do you want to Solve This Problem? (1=yes, 0=no). 0
syntax highlighted by Code2HTML, v. 0.8.8b
with associated plot
