Initial Value Problems - Numerical Methods in Engineering with MATLAB - page 264

Solution The function specifying the differentialequationis

_

functionF=fex7 5(x,y)

% Differential eq. used in Example 7.5.
F=3*y-4*exp(-x);

The solutionis(every 20th line was printed):

>> [x,y] = runKut4(@fex7

_

5,0,1,10,0.1);

>> printSol(x,y,20)
x y1
0.0000e+000 1.0000e+000
2.0000e+000 1.3250e-001
4.0000e+000 -1.1237e+000
6.0000e+000 -4.6056e+002
8.0000e+000 -1.8575e+005
1.0000e+001 -7.4912e+007

It is clear thatsomething went wrong. According to the analytical solution, y
shoulddecrease to zero with increasing x , but the output shows the opposite trend:
after an initial decrease, the magnitudeof y increases dramatically. The explanation
isfound by taking a closer look at the analytical solution. The general solution of the
givendifferentialequationis

Ce 3 x e − x

= +

y

which can be verifiedbysubstitution. The initialcondition y (0) = 1yields C = 0, so

e − x .
The cause of trouble in the numerical solutionis the dormantterm Ce 3 x . Suppose
that the initialcondition contains a small error

that the solution to the problemis indeed y =

ε ,sothat wehave y (0) = 1 + ε . This

changes the analytical solution to

e 3 x e − x y = ε +

We now see that the term containing the error becomes dominant as x is increased.

Since errors inherent inthenumerical solutionhave the same effect assmall changes in
initialconditions, weconcludethatour numerical solutionis the victim of numerical
instability duetosensitivity of the solution to initialconditions. The lesson here is:do
not alwaystrust the results of numerical integration.

ε