MATLAB Answers

How to find unstable solution

5 views (last 30 days)
Betty Johnson
Betty Johnson on 8 Oct 2020
Commented: Betty Johnson on 8 Oct 2020
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.

  2 Comments

Mathieu NOE
Mathieu NOE on 8 Oct 2020
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
Betty Johnson on 8 Oct 2020
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.

Sign in to comment.

Answers (1)

John D'Errico
John D'Errico on 8 Oct 2020
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 = 
et
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 = 
et(75e12t+9007199254741817)9007199254740992
vpa(expand(ysol2),20)
ans = 
1.0000000000000915934e1.0t+0.0000000000000083266726846886740532e11.0t
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.

  1 Comment

Betty Johnson
Betty Johnson on 8 Oct 2020
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.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by