My Favorite ODE
My favorite ordinary differential equation provides a good test of ODE solvers, both numeric and symbolic. It also provides a nice illustration of the underlying existence theory and error analysis.
Contents
Implicit equation
Here is one of my favorite ordinary differential equations. The dot means differentiation with respect to t. The equation does not express ˙y directly, so it is implicit.
˙y2+y2=1
We don't need a computer to verify that the analytic solution is a combination of sint and cost. But the constants y=1 and y=−1 are also solutions.
Explicit equation
Our ODE solvers require an explicit equation, but when we solve this equation for ˙y we get two solutions, differing in sign.
˙y=±√1−y2
That ± sign makes this interesting.
First attempt
Let's plunge ahead with the + sign. Specify the equation, an interval for t and an initial condition.
f = @(t,y) sqrt(1-y^2)
format short
tspan = [0 pi]
y0 = 0
f = @(t,y)sqrt(1-y^2) tspan = 0 3.1416 y0 = 0
An attempt to compute a numeric solution with ode45 runs into trouble when y gets slightly larger than 1 and the square root of a negative number generates an imaginary result. So temporarily use a shorter t interval.
try ode45(f,tspan,y0); catch ode45(f,[0 1.47],y0) end pi_axis([-0.1 1.1])
Error using odeplot (line 63) Error updating the ODEPLOT window. Solution data may have been corrupted. Argument Y cannot be complex. Error in ode45 (line 435)
Keep in bounds
Change the equation so that it provides a guard against the square root of negative numbers.
f = @(t,y) sqrt(1-min(y,1)^2) ode45(f,tspan,y0); pi_axis([-0.1 1.1])
f = @(t,y)sqrt(1-min(y,1)^2)
This is the correct solution, sint. For a while. But when y reaches 1 it has to stay there because its derivative is always nonnegative.
Switch sign
If we want to continue with sint beyond π/2 we have to switch to the minus sign.
f = @(t,y) (1-2*double(t>pi/2))*real(sqrt(1-y.^2)) ode45(f,tspan,y0) pi_axis([-0.1 1.1])
f = @(t,y)(1-2*double(t>pi/2))*real(sqrt(1-y.^2))
Jacobian
We're solving the autonomous differential equation,
˙y=f(y)
where
f(y)=√1−y2
An important quantity in both the existence theory for ODEs and the error analysis for numeric methods is the Jacobian.
J=∂f∂y
In this case J is 1-by-1.
J=−y√1−y2
When y approaches ±1 the Jacobian approaches infinity.
Existence and uniqueness theorems for ODEs rely on a bounded Jacobian (or a bounded Lipschitz constant). The derivations of numerical methods rely on a bounded Jacobian. Error control in modern ODE software relies on a bounded Jacobian. All of these are up in the air with this example.
Nonuniqueness
We can see the lack of uniqueness happening when we change the initial condition to y(0)=−1. The Jacobian is infinite at the start and both the constant y(t)=−1 and the trig function y(t)=−cost are solutions. A small perturbation in the initial condition switches the numerical solution from one to the other.
f = @(t,y) sqrt(1-max(min(y,1),-1)^2) y0 = -1 ode45(f,tspan,y0); pi_axis([-1.2 0]) hold on format long y0 = -.99999999 ode45(f,tspan,y0); pi_axis([-1.1 1.1]) hold off
f = @(t,y)sqrt(1-max(min(y,1),-1)^2) y0 = -1 y0 = -0.999999990000000
Symbolic solutions
I was curious to see how the Symbolic Toolbox handles this equation with this initial condition. There is good news and surprising news.
syms y(t)
ode = diff(y)^2+y^2==1
init = y(0)==-1
ysym = dsolve(ode,init)
ode(t) = diff(y(t), t)^2 + y(t)^2 == 1 init = y(0) == -1 ysym = -1 -cosh(t*1i)
I am pleased to see the two solutions confirmed, but I am surprised to see the second one expressed as a hyperbolic trig function of an imaginary argument.
Simplify
Simplification generates a more familiar form.
simplify(ysym)
ans = -1 -cos(t)
Second order equation
Let's return to the original implicit equation and differentiate the entire equation with respect to t.
ddt(˙y2+y2=1)
2˙y¨y+2y˙y=0
We can divide through this equation by 2˙y, provided that ˙y is nonzero. This eliminates the constant solutions and produces an old friend, the harmonic oscillator.
¨y+y=0
First order system.
Our ODE solvers need the second order equation to be expressed as a first order system.
˙y1=y2, ˙y2=−y1
Now there is no trouble computing sint and −cost.
f = @(t,y)[y(2); -y(1)] y0 = [0; -1] ode23(f,tspan,y0) pi_axis([-1 1])
f = @(t,y)[y(2);-y(1)] y0 = 0 -1
pi_axis
I have been using this little function to label the x-axis with multiples of π and scale the y-axis.
type pi_axis
function pi_axis(yaxis) axis([0 pi yaxis]) set(gca,'xtick',0:pi/4:pi, ... 'xticklabels',{0 '\pi/4' '\pi/2' '3\pi/4' '\pi'}) end
댓글
댓글을 남기려면 링크 를 클릭하여 MathWorks 계정에 로그인하거나 계정을 새로 만드십시오.