ORDINARY DIFFERENTIAL EQUATIONS POWERTOOL

Unit 23 -- Nonlinear Equations

Outline of Unit 23

Initialization

> restart;

> with( DEtools ):

> with( plots ):

> with( linalg ):

> with( student ):

> with( plottools ):

Warning, the name changecoords has been redefined

Warning, the name adjoint has been redefined

Warning, the protected names norm and trace have been redefined and unprotected

Warning, the name translate has been redefined

> COLORS := [ BLUE, GREEN, CYAN, AQUAMARINE ]:

>

To solve a nonlinear second-order ODE of the form

> ODE := diff( y(x), x$2 ) = F( x, diff( y(x), x ) );

>

introduce a second dependent function

> new_dep_var := w(x) = diff( y(x), x );

>

Then = = and the ODE can be rewritten as a first-order system

> sys_ode1 := isolate( new_dep_var, diff( y(x), x ) );

> sys_ode2 := eval( ODE, sys_ode1 );

>

The second-order ODE for has become a first-order ODE for . Assuming this ODE can be solved (for ) by one of the first-order methods discussed previously, integration of the first-order ODEf for , i.e. , , provides a solution to the second-order ODE.

>

Consider the second-order nonlinear ODE

> ode1 := 2*x*diff(y(x),x$2) - diff(y(x),x) + 1/diff(y(x),x) = 0;

>

for > 0 with "initial conditions"

> ic1 := y(1) = 2, D(y)(1) = -2;

>

Since this equation does not contain the dependent variable, , re-write it as the first-order system

> sys1 := [ sys_ode1, eval( ode1, sys_ode1 ) ];

>

with initial conditions

> sys_ic1 := eval( ic1, D(y)=w );

>

Maple classifies the ODE for as a separable equation.

> odeadvisor( sys1[2] );

>

and finds the solution to be

> sol_w := dsolve( { sys1[2], sys_ic1[2] }, w(x), [separable] );

>

Now, the remaining first-order ODE is

> q1 := eval( sys1[1], sol_w );

>

which can be solved by direct integration

> sol1 := dsolve( { q1, sys_ic1[1] }, y(x), [quadrature] );

>

Checking that this function does solve the second-order ODE

> odetest( sol1, ode1 );

>

and initial conditions

> simplify( eval( sol1, x=1 ) );

> simplify( eval( q1, x=1 ) );

>

Also, Maple reports the same solution when dsolve is used to solve the second-order IVP

> infolevel[dsolve] := 3:

> dsolve( { ode1, ic1 }, y(x) );

> infolevel[dsolve] := 0:

Methods for second order ODEs:

Trying to isolate the derivative d^2y/dx^2...

Successful isolation of d^2y/dx^2

-> Trying classification methods

trying 2nd order Liouville

trying 2nd order, 2 integrating factors of the form mu(x,y)

trying differential order: 2; missing variables

Found another symmetry   [x, y]

   *** Sublevel 2 ***

   symmetry methods on request

   1st order, trying reduction of order with given symmetries:   [_a, 0]

   1st order, trying the canonical coordinates of the invariance group

   1st order, canonical coordinates successful

differential order: 2; canonical coordinates successful

differential order 2; missing variables successful

>

The solution procedure for a second-order ODE that does not depend explicitly on the independent variable

> ODE2 := diff( y(x), x$2 ) = f( y(x), diff( y(x), x ) );

>

begins with the introduction of the new dependent variable

> new_dep_var;

>

which allows the second-order ODE to be written as the two equations

> sys_ode1;

> q2 := eval( ODE2, sys_ode1 );

>

This form is not yet useful as the independent variable is still involved, albeit implicitly, in these equations. To complete the transformation to a useful form, change the dependent variable from x to y. That is, use the Chain Rule to write

= = .

> new_indep_var2 := diff( w(x), x ) = diff( w(y), y ) * w(y);

> new_dep_var2 := y(x)=y, w(x)=w(y);

>

The ODE for with independent variable can now be written as

> ODE2_w := eval( q2, { new_indep_var2, new_dep_var2 } );

>

Assuming this first-order ODE can be solved explicitly, the solution to the first-order ODE is

> ODE2_y := eval( sys_ode1, { new_dep_var2[2] } );

>

is the solution to the original second-order ODE.

Consider the nonlinear second-order ODE

> ode2 := diff( y(x), x$2 ) - diff( y(x), x )^3 = 0;

>

The changes of dependent and independent variable is completed in two steps

> ode2_wx := eval( ode2, sys_ode1 );

> ode2_wy := eval( ode2_wx, { new_indep_var2, new_dep_var2 } );

>

The separable ODE for has general solution

> sol2_wy := dsolve( ode2_wy, w(y), [separable] );

>

To complete the solution process, invert the change of independent variable

> sol2_wx := subs( [ w(y)=w(x), y=y(x) ], sol2_wy );

>

to obtain the ODE for

> ode2_yx := eval( sys_ode1, sol2_wx );

>

The general (implicit) solution to this separable ODE is

> sol2 := dsolve( ode2_yx, y(x), [separable], implicit );

>

Initial conditions are required to determine the constants _C1 and _C2 . In lieu of these data values, it is easily checked that these solutions satisfy the second-order ODE

> odetest( sol2, ode2 );

>

Recall that the qualitative analysis of solutions to a nonlinear ODE is easiest to discuss after the ODE has been converted to a first-order system. Then, the behavior of solutions near an equilibrium point is determined by the sign of the eigenvalues of the Jacobian matrix of the system at the equilibrium point.

> X := [ x(t), y(t) ]:

> U := [ u(t), v(t) ]:

>

Consider the nonlinear system

> f := (x,y) -> -3*x + 3*x^2:

> g := (x,y) -> -2*x + 2*x^2 + 3*y - 3:

> sys := equate( diff( X, t ), [ f(x(t),y(t)), g(x(t),y(t)) ] );

>

The direction field for this nonlinear system

> Pdfield := DEplot( sys, [ x(t), y(t) ], t=0..1, x=-2..2, y=-2..2,

> arrows=SLIM, color=PINK ):

> display( Pdfield );

>

suggests that there are (at least) two equilibrium points for this system

> equil := solve( {f(x0,y0)=0,g(x0,y0)=0}, {x0,y0} );

>

> Pequil := pointplot( map( subs, [equil], [x0,y0] ),

> symbol=CIRCLE, symbolsize=18, color=RED ):

> display( Pdfield, Pequil );

>

Linearization of the system about each of these equilibrium points

> for i from 1 to nops([equil]) do

> X0 := equil[i];

> b := eval( [ f(x0,y0), g(x0,y0) ], X0 );

> A := eval( evalm(jacobian( [f(x0,y0),g(x0,y0)], [x0,y0] )), X0 );

> linsys1 := equate( diff( U, t ), evalm( A &* U + b ) );

> Plin||i := translate( DEplot( linsys1, {u(t),v(t)}, t=0..1,

> u=-1/2..1/2, v=-1/2..1/2,

> arrows=SLIM, color=COLORS[i],

> dirgrid=[7,7] ),

> op(eval([x0,y0],X0)) );

> print( `-------------------------------------------------------------` );

> print( `Linearization at [x0,y0]`=eval([x0,y0],X0) );

> print( Plin||i );

> print( nprintf( "%s %a", `Eigenvalues:`, [eigenvals( A )] ) );

> print( `-------------------------------------------------------------` );

> end do:

>

Thus, the equilibrium solution is a saddle point and is an unstable node.

To conclude, superimpose the direction fields for the linearizations on the original direction field.

> display( Pdfield, Pequil, Plin||(1..nops([equil])), view=[-1..2, 0..2 ] );

>

Consider the nonlinear system

> f := (x,y) -> (2-x)*x + x*y/2:

> g := (x,y) -> (3-y)*y + x*y/2:

> sys := equate( diff( X, t ), [ f(x(t),y(t)), g(x(t),y(t)) ] );

>

The direction field for this nonlinear system

> Pdfield := DEplot( sys, [ x(t), y(t) ], t=0..1, x=-1..7, y=-1..7,

> arrows=SLIM, color=PINK ):

> display( Pdfield );

>

suggests that there are (at least) four equilibrium points for this system

> equil := solve( {f(x0,y0)=0,g(x0,y0)=0}, {x0,y0} );

>

> Pequil := pointplot( map( subs, [equil], [x0,y0] ),

> symbol=CIRCLE, symbolsize=18, color=RED ):

> display( Pdfield, Pequil );

>

Linearization of the system about each of these equilibrium points

> for i from 1 to nops([equil]) do

> X0 := equil[i];

> b := eval( [ f(x0,y0), g(x0,y0) ], X0 );

> A := eval( evalm(jacobian( [f(x0,y0),g(x0,y0)], [x0,y0] )), X0 );

> linsys1 := equate( diff( U, t ), evalm( A &* U + b ) );

> Plin||i := translate( DEplot( linsys1, {u(t),v(t)}, t=0..1,

> u=-1/2..1/2, v=-1/2..1/2,

> arrows=SLIM, color=COLORS[i],

> dirgrid=[7,7] ),

> op(eval([x0,y0],X0)) );

> print( `-------------------------------------------------------------` );

> print( `Linearization at [x0,y0]`=eval([x0,y0],X0) );

> print( Plin||i );

> print( nprintf( "%s %a", `Eigenvalues:`, [evalf(eigenvals( A ))] ) );

> print( `-------------------------------------------------------------` );

> end do:

>

Thus, the equilibrium solution is an unstable node (source), is a saddle point, is a saddle point, and is a stable node (sink).

To conclude, superimpose the direction fields for the linearizations on the original direction field.

> display( Pdfield, Pequil, Plin||(1..nops([equil])), view=[-1..7, -1..7 ] );

>