Names:

Load these packages:

> map(with,{plots,DEtools}):

Section 2.3 of the Borrelli and Coleman textbook discusses the sensitivity of solutions to IVPs

under changes to the data, f and a. Roughly speaking, do small changes to the driving function, f, and or the initial data point, a, produce small (or even predictable) changes to the solutions y(t)?

For starters, we simply "plug in different values" and see what happens. Consider the IVP

and study the changes induced by altering slightly the value of the parameter c. Using the trick of treating c as a state variable (c(t) = constant, or c'(t) = 0) in the text, we can enter this IVP as a system:

> ode1:={diff(y(t),t)=-y(t)*sin(t) + c(t)*t*cos(t), diff(c(t),t) = 0};

The initial condition c(-4) = 5, guarantees or forces c(t) = 5 for all t, so to "trick" Maple into plotting solutions for various values of c, we simply feed DEplotseveral initial conditions. Generate the list of initial conditions (from the text):

> inits:=[seq([y(-4)=-6, c(-4)=.7+k*.1], k=0..6)];

This is essentially the second plot on the top of page 108.

> DEplot(ode1,[y(t),c(t)], t=-4..4, inits, y=-10..15,
scene=[t,y],linecolor=black);

We see that small changes in c produce large variations in the solution near t = 0, but then (perhaps) the solutions come back together again for larger values of t. Near t = -4 (the initial point), all the solutions appear to be close together. The bottom solution corresponds to c = 0.7 and the top solution corresponds to c = 1.3. Let's see if the curves for c = 0 .7 and c = 1.3 "run together" as t increases beyond t = 4. I'll use the dsolve(..., numeric, ...) techniques from the lab-book, chapter 7, to create procedures to give values for these solutions, and then plot the solutions. Since we are dealing with procedures, the plot range is specified without the var =part of the usual plot command.

> ansLow:=dsolve({diff(y(t),t)=-y(t)*sin(t) + 0.7*t*cos(t),y(-4)=-6},
y(t), numeric,startinit=true):

ansHigh:= dsolve({diff(y(t),t)=-y(t)*sin(t) + 1.3*t*cos(t),y(-4)=-6},
y(t), numeric, startinit=true):

yLow:= u -> subs(ansLow(u),y(t)):

yHigh:= u-> subs(ansHigh(u),y(t)):

> plot([yLow, yHigh],-4..10,color=[red,blue]);

The numeric approximations begin to show the long term effects of changing the value of c. Except near t = 0, the solutions (to different IVPs, since the value of c is different) have a lot in common, but are somewhat different. Typically, near the initial point t = -4, solutions should look highly similar, far away from t = -4, there is no reason to expect similarities in solutions to slightly different IVPs. The fact that these two different solutions are so similar from t = 4 to 8 is abnormal. Here are "formulas" for the bottom (red) and top (blue) solutions plotted above:

> dsolve({diff(y(t),t)=-y(t)*sin(t) + 0.7*t*cos(t),y(-4)=-6},
y(t));

dsolve({diff(y(t),t)=-y(t)*sin(t) + 1.3*t*cos(t),y(-4)=-6},
y(t));

You can see that the constants, c = 0.7 and c = 1.3, are multiplied by a fairly complicated expression (an exponential of cosine times an area function), as the expression which is multiplied by c grows, so does the difference between the 2 solutions. Here is a plot of the expression by which c is multiplied:

> plot(exp(cos(t))*int(u*cos(u)*exp(-1.*cos(u)),u = -4 .. t),
t=-5..15);

Near t = -4, the multiplier is small, and thus the solutions will be nearly the same. Likewise, for values of t from 3 to 6, or near t = 8.5, the multiplier is small. On these intervals, we say that the IVP is insensitive to variations in the parameter c.

From t = -2 to 2 and for t > 9, the multiplier gets quite large, and the 2 solutions will differ significantly. This means that the IVP is sensitive to variations in the parameter c on these intervals.

We next look at Theorem 2.3.1, which involves the linear IVP y' + p(t)y = q(t), y(0) = a, where p(t) and q(t) are continuous on the interval in question, say I = all t from 0 to 10. The assumptions state that there are positive constants P and M with p(t) >= P and | q(t) | <= M (for all t in I). This is known as Bounded Inputor BI. The conclusions state that any (solutions are unique here) solution is trapped between the constants | a | + M / P and -| a | - M / P, for t > 0. Likewise, a tighter bound on the absolute value of a solution is given by

where P can be negative in this case.

Consider y' + ((9 + t) / (3 + t))*y = 5 - 3*exp(-t), y(0) = 2. Here are plots of p(t) (red) and q(t) (blue):

> p:= t -> (9+t)/(3+t):
q := t -> 5 - 3*exp(-t):

> plot([p(t),q(t)],t=0..10,y=0..5, color=[red,blue]);

In order to use the theorem's bounds, we need to find a lower bound, say P, for p(t), where P <= p(t) for all t in our interval. The global minimum of p(t) on the interval I is the best choice for P. The plot strongly suggests that p(10) is the global minimum.

> p(10.);

Any value of P smaller than or equal to that value will do. I'll use P = 1.46, as it makes sense to take P close to the global minimum.

We also need an upperbound, say M, for q(t), with q(t) <= M for all t in the interval I. The global maximum of q(t) on I is a good choice for M, and this appears to occur at t = 10 as well. Recall that global extrema occur either at endpoints of the interval in question, or at values of t where the derivative is zero or undefined, i.e., p'(t) = 0 or q'(t) = 0.

> q(10.);

Any value bigger than or equal to that will do for M, I'll use M = 5.

Using P = 1.46, M = 5 and a = 2 in the theorem, gives the following functional and constant bounds.

> highf:= t-> exp(-1.46*t)*2 + (5/1.46)*(1-exp(-1.45*t));
lowf:= t -> -highf(t);
highc:= 2+5/1.46;
lowc:=-highc;

A plot which defines and then shows these "bounds" as thick black curves:

> bnds:=plot([highf(t), lowf(t), highc, lowc],
t=0..10,color=black,thickness=3):

bnds;

Here is the numeric solution along with these bounds:

> display([bnds, DEplot(diff(y(t),t)+p(t)*y(t)=q(t),
y(t),t=0..10, [[y(0)=2]], y=-8..8, linecolor=blue,
stepsize=.05)]);

It was predicted in advance that the blue curve would lie between the bounding constants ( y = +- (a + M/P)) and the bounding functions. The power here is that one only requires knowledge that p(t) is bigger than some fixed value (P), that q(t) is smaller than some fixed value (M), and one can then control the size of the solution based solely on these bounds. Using the global minimum of p(t) and the global maximum of q(t) (on the interval I) is a good choice.

One more example with much less explanation. Consider y' + (3 + t*ln(t + 0.2) )*y = 1 + t*cos(t), y(0) = 1, for 0 <= t <= 2. The main idea is to define p(t) and q(t), locate their local extrema, calculate P and M, then plot it all. I start by "undefining" p and q, to avoid confusing the old values above, with the newer values here.

> p:= 'p':
q:= 'q':

p:= t -> 3+t*ln(t+.2);
q:= t -> 1 + t*cos(t);

> plot([p(t),q(t)],t=0..2,color=[red,blue]);

Global extrema occur where the derivative is zero. The plot provides my "guesses".

> tp:= fsolve(diff(p(t),t) = 0, t=0.3);
tq:= fsolve(diff(q(t),t) = 0, t=0.9);

> P := p(tp);
M := q(tq);
a := 1.;

> highf:= t-> exp(-P*t)*a + (M/P)*(1-exp(-P*t));
lowf:= t -> -highf(t);
highc:= a+M/P;
lowc:=-highc;

> bnds2:=plot([highf(t), lowf(t), highc, lowc],
t=0..2,color=black,thickness=3):

bnds2;

> display([bnds2, DEplot(diff(y(t),t)+p(t)*y(t)=q(t),
y(t),t=0..2, [[y(0)=a]], y=-1.8..1.8, linecolor=blue,
stepsize=.05)]);

Execute this command to clear the definitions on P, M, a, p and q.

> map(unassign,{'P','M','a','p','q'}):

Justify the claim that p(t) = 2 + cos(t) and q(t) = 4.8*sin(t) satisfy input bounds (P = ? and M = ?) similar to the problems above, on the interval t = 0 to 10, then make a plot of the solution to this new IVP (and the new bounds), with y(0) = 3, similar to the plots above.
 
 

Theorem 2.3.2 addresses situations where the initial value a = y(0) and the input function, q(t), for the IVP

are subject to small (often uncontrolable) variations. If we call the actual solution z(t), the actual initial value b (so in practice, z(0) = b is the true initial value) and the actual input function m(t) (so the actual solution satisfies z' + p(t)*z = m(t), z(0) = b), the hope is that we can control discrepencies between the ideal solution y(t) of (*) and the actual solution z(t). Under the assumptions that p(t), q(t) and m(t) are all continuous, and the existence of positive constants P and M, where

the deviation between y(t) and z(t) is bounded by

Consider the ideal IVP y' + 0.6y = 3.5*sin(2*t), y(0) = 1. Here is the ideal solution

> ysoln:=rhs(
dsolve({diff(y(t),t)+0.6*y(t) = 3.5*sin(2*t), y(0)=1},y(t)));

If in actuality, we know the initial value is between 0.9 and 1.1 (so | a - b | <= 0.1), and we add a small "wiggle" to q(t), setting m(t) = q(t) + 0.2*cos(t), then | q(t) - m(t) | <= | 0.2*cos(t) | <= 0.2, since | cos(t) | <= 1, we can take P = 0.6, and M = 0.2 in the theorem. Here is a plot of the maximum possible deviations:

> dev:=exp(-0.6*t)*0.1 + (0.2/0.6)*(1-exp(-0.6*t));

> plot([dev,-dev],t=0..10,color=red);

Now, add these maximal deviations to the ideal solution:

> band:=
plot([ysoln+dev,ysoln-dev],t=0..10,color=red):

band;

All solutions satisfying the bounding conditions on m(t) and b, must lie inside that band. Here is the ideal solution along with the band of maximal error:

> plot([ysoln,ysoln+dev,ysoln-dev],t=0..10,color=[blue,red,red]);

We see in the plot above that any actual solution remains close to the ideal solution for all t <= 10 (and larger t values actually stay close in this case also). Here is the solution (numeric) corresponding to b = 0.92 and the wiggled version m(t) = q(t) + 0.2cos(t), from above.

> zans:=dsolve(
{diff(z(t),t)+0.6*z(t) = 3.5*sin(2*t)+0.2*cos(t),
z(0)=0.95},z(t),numeric):

zsoln:= u-> subs(zans(u),z(t)):

And a graph of this actual solution along with the band:

> display([plot(zsoln,0..10,color=black), band]);

Assume instead that | m(t) - q(t) | <= 0.3*t, for 0 <= t <= 4, and that z(0) = b can differ from 1 by as much as 0.25 (in absolute value). Produce an error band plot (centered at the ideal solution named ysoln above), and then graph the solution, inside the band, which satisfies z' + 0.6z = 3.5*sin(2*t) + 0.29*t, z(0) = 1.23.
 
 

Putting it all together:

The solution of y' + y = c*t, y(0) = a is denoted by y_c(t, a). We want to look at changes in these solutions as we vary the values of c and a. Problem 5c asks us to investigate these questions in comparison to c = 1, for t between 0 and 10. This IVP is linear with p(t) = 1, and q(t) = c*t. Since we are allowed only to change the values of a and c, we are forced to consider input functions m(t) = d*t, and initial conditions z(0) = b.

For starters, look at a DEplot, with c = 1 and a variety of values for b, say 5 values close to a = 4.
 
 

For t from 0 to 10, does it appear that this IVP, with b near 4, is sensitive or insensitive to variations in the initial point b?

Using the trick of making c a state variable (add c'(t) = 0 to the system, etc.), make a plot of the solutions with z(0) = 4 and c = 0.5, 0.6, ... , 1.5.
 
 

Does it appear that this IVP (with a = 4) is sensitive or insensitive to variations in the input function q(t)?
 
 

By solving the differential equation, find a formula for (t,a).
 
 

Show (directly, without using theorem 2.3.2) that

and use this inequality to explain why (t, a) is not sensitive to variations in a, but is highly sensitive to variations in c, for t from 0 to 10.
 
 

Find good values for P and M (where z(t) = (t,b) ), and use theorem 2.3.2 to produce a bound, similar to the one above, for | (t,a) - (t,b) |.