How to find unstable solution

How to get two lines in one graph for

y''-10y'-11y=0 with initial conditions y(0)=1,y'(0)=-1

and with the change in the first intial condition y(0)=1+epsilon ,y'(0) =-1

I run Heun's method but the solutions are slightly different though it should show the stability for the first line and unstability for the 2nd line.

[Mathieu NOE]

hello betty

IMHO, your second order ode is always unstable , whatever the IC

are you sure about the sign of the coefficients ?

demo : let Y = [y' ; y]

your 2nd order ode can be rewritten as a state space representation :

Y' = A Y

with A = [10 11; 1 0]

you compute the eigen values of A ; there is one eigenvalue positive so system is unstable (to be stable , all eigenvalues must be strictly negative)

[v,d] = eig(A)

d =

11 0

0 -1

in my simulations, I always got diverging curves

you may be fooled if you plot both curves in the same plot,

because one diverges much faster than the other one (so due to y axis scale, you'd thought that the less diverging curve is stable)

[Betty Johnson]

Thank you for your response.I used the following codes and I get almost same answer as yours but it seems like that the graph should look like the attached figure.

Answer 1

[John D'Errico]

Your question about how to get two lines in one graph is a simple one. Use the command "hold on" between plots. That is...

plot(1:5,rand(1,5),'r-o')
hold on
plot(1:5,rand(1,5),'b-x')
hold off

So now we have one plot, with two curves drawn on it.

Now look at the problem you are trying to solve. Though I won't do your homework for you, writing a Heun's method solver to produce several results plotted together, we can look at the problem and explain what should happen. It was:

y''-10y'-11y=0 with initial conditions y(0)=1,y'(0)=-1

What is happening there? Look at this curve, in the form of a symbolic solution.

syms y(t)
dy = diff(y,t);
ysol = dsolve(diff(y,t,2) - 10*dy - 11*y == 0,y(0) == 1,dy(0) == -1)

ysol = 

e−t

As you should see, this solution is indeed a nice, simple thing. In the analytical solution, we have a second term in there multiplying exp(11*t), but in the analytical solution for the unperturbed problem the multipler of exp(11*t) was exactly ZERO. The full analytical solution could have been written

ysol = exp(-t) + 0*exp(11*t)

given that set of initial conditions. That is theory. Theory is always so nice, right? But then practice gets in the way.

What happens now if we change the initial conditions just by a tiny amount? A tiny error in the nether-bits is enough to nudge the solver to give me a second non-zero term in there.

ysol2 = dsolve(diff(y,t,2) - 10*dy - 11*y == 0,y(0) == 1.0000000000001,dy(0) == -1)

ysol2 = 

e−t (75 e12 t+9007199254741817)9007199254740992

vpa(expand(ysol2),20)

ans = 

1.0000000000000915934 e−1.0 t+0.0000000000000083266726846886740532 e11.0 t

So what happens with this set of initial conditions is as the solver moves along, you need to remember you are using a numerical ODE solver. Is the solution produced by any numerical solver exact? No. In fact, you should know it has an error term. You should expect some tiny amount of slop in the solution. But now once we wander off the true colution, even by the tiniest amount, that second part comes into the problem, and it introduces itself with a vengeance. That second term gets HUGE, and it does so very rapidly. So as soon as any tiny amount of crap gets into the problem, the ODE solver can start to chase along the bad solution.

It will get worse if you start away from the stable part of the solution by just a tiny amount. Now the solver is already looking at something that will get huge quickly. It does not even need to introduce numerical errors into the problem, to then be massively amplified in later steps.

[Betty Johnson]

Hi,thank you for your response.

I used the following codes but it doesn't give me any different graph by changing the initial condition,though I defined it a way that it should change.