rifsimp - simplify overdetermined polynomially nonlinear PDE or ODE systems

Calling Sequences

rifsimp( system , options )

rifsimp( system , vars , options )

Parameters

system - list or set of polynomially nonlinear PDEs or ODEs (may contain inequations)

vars - (optional) list of the main dependent variables

options - (optional) sequence of options to control the behavior of rifsimp

Description

Examples

1. Overdetermined Systems

As a first example, we have the overdetermined system of two equations in one dependent variable y(x) .

> sys1:=[x*diff(y(x),x)^2*diff(y(x),x,x)^2-2*x*diff(y(x),x)*
diff(y(x),x,x)*y(x)*diff(y(x),x,x,x)+x*y(x)^2*
diff(y(x),x,x,x)^2-y(x)*diff(y(x),x,x)+diff(y(x),x)^2=0,
-diff(y(x),x)*diff(y(x),x,x)+y(x)*diff(y(x),x,x,x)+
2*y(x)^2*diff(y(x),x,x)^2-4*y(x)*diff(y(x),x,x)*
diff(y(x),x)^2+2*diff(y(x),x)^4=0];

sys1 := [x*diff(y(x),x)^2*diff(y(x),`$`(x,2))^2-2*x...
sys1 := [x*diff(y(x),x)^2*diff(y(x),`$`(x,2))^2-2*x...
sys1 := [x*diff(y(x),x)^2*diff(y(x),`$`(x,2))^2-2*x...
sys1 := [x*diff(y(x),x)^2*diff(y(x),`$`(x,2))^2-2*x...
sys1 := [x*diff(y(x),x)^2*diff(y(x),`$`(x,2))^2-2*x...

Now call rifsimp .

> with(Rif):
ans1:=rifsimp(sys1);

ans1 := TABLE([Solved = [diff(y(x),`$`(x,2)) = diff...
ans1 := TABLE([Solved = [diff(y(x),`$`(x,2)) = diff...
ans1 := TABLE([Solved = [diff(y(x),`$`(x,2)) = diff...
ans1 := TABLE([Solved = [diff(y(x),`$`(x,2)) = diff...

We see that the system has been reduced (for the case y(x) not identically zero) to one much simpler equation that can now be handled with Maple's dsolve .

> dsolve(convert(ans1[Solved],'set'),{y(x)});

{y(x) = exp(_C1*x)*_C2}

In addition to assisting dsolve in obtaining exact solutions, the simplified form can be used in combination with the initialdata function to obtain the required initial data for the solved form.

> id1:=initialdata(ans1[Solved]);

id1 := TABLE([Finite = [y(x[0]) = _C1, diff(y(x[0])...

In this case, it gives us the expected results for a second order ODE, but it can also calculate the required initial data for more complex PDE systems (see initialdata for more information). Numerical methods can now be successfully applied to the reduced system, with initial conditions of the type calculated by initialdata .

> dsolve(convert(ans1[Solved],set) union {y(0)=1,D(y)(0)=2},
{y(x)}, numeric, output=array([0,0.25,0.5,0.75,1.0]));

Error, (in dsolve/numeric/an_args/common) 'output' argument
must be one of, 'procedurelist','listprocedure'

Finally, we can also obtain a formal power series (Taylor series) for the local solution of this problem (see rtaylor for more information).

> tay_ser:=rtaylor(ans1[Solved],order=4);

tay_ser := [y(x) = y(x[0])+diff(y(x[0]),x)*(x-x[0])...
tay_ser := [y(x) = y(x[0])+diff(y(x[0]),x)*(x-x[0])...
tay_ser := [y(x) = y(x[0])+diff(y(x[0]),x)*(x-x[0])...

The above series is fully determined when given the initial data described earlier:

> tay_ser:=rtaylor(ans1[Solved],id1,order=4);

tay_ser := [y(x) = _C1+_C2*(x-x[0])+1/2*_C2^2/_C1*(...
tay_ser := [y(x) = _C1+_C2*(x-x[0])+1/2*_C2^2/_C1*(...

2. Inconsistent Systems

rifsimp can often determine whether a system is inconsistent.

> sys2:=[diff(u(x,y),x)^2+diff(u(x,y),y,y)=0,
u(x,y)^3-diff(u(x,y),x)+x=0]:
rifsimp(sys2);

TABLE([status =

3. Constrained Mechanical Systems

This example shows the use of rifsimp as a preprocessor for a constrained mechanical system (that is, a Differential-Algebraic Equation or DAE system). The method of Lagrange formulates the motion of a bead of nonzero mass m on a frictionless wire of shape Phi(x,y) under gravity as follows:

> ConstrainedSys:= [m<>0, Phi(x(t),y(t))=0,
m*diff(x(t),t,t)=-lambda(t)*D[1](Phi)(x(t),y(t)),
m*diff(y(t),t,t)=-lambda(t)*D[2](Phi)(x(t),y(t))-m*g];

ConstrainedSys := [m <> 0, Phi(x(t),y(t)) = 0, m*di...
ConstrainedSys := [m <> 0, Phi(x(t),y(t)) = 0, m*di...
ConstrainedSys := [m <> 0, Phi(x(t),y(t)) = 0, m*di...

Note that a mass falling under gravity without air resistance corresponds to:

> Phi:= (x,y) -> 0;

Phi := 0

We pick the example of a pendulum, so the bead moves on a circular wire.

> Phi:= (x,y) -> x^2 + y^2 - 1;

Phi := proc (x, y) options operator, arrow; x^2+y^2...

> ConstrainedSys;

[m <> 0, x(t)^2+y(t)^2-1 = 0, m*diff(x(t),`$`(t,2))...
[m <> 0, x(t)^2+y(t)^2-1 = 0, m*diff(x(t),`$`(t,2))...

Such constrained systems present great difficulties for numerical solvers. Of course, we could eliminate the constraint using polar coordinates, but in general this is impossible for mechanical systems with complicated constraints. For example, just replace the constraint with another function, so that the wire has a different shape.

> Phi2:= (x,y) -> x^4 + y^4 - 1;

Phi2 := proc (x, y) options operator, arrow; x^4+y^...

(See Reid et al. , 1996, given in the reference and package information page Rif ). The pendulum is the classic example of such systems that have been the focus of much recent research due to their importance in applications.

> SimpSys:=rifsimp(ConstrainedSys);

SimpSys := TABLE([Solved = [diff(y(t),`$`(t,2)) = -...
SimpSys := TABLE([Solved = [diff(y(t),`$`(t,2)) = -...
SimpSys := TABLE([Solved = [diff(y(t),`$`(t,2)) = -...
SimpSys := TABLE([Solved = [diff(y(t),`$`(t,2)) = -...
SimpSys := TABLE([Solved = [diff(y(t),`$`(t,2)) = -...
SimpSys := TABLE([Solved = [diff(y(t),`$`(t,2)) = -...
SimpSys := TABLE([Solved = [diff(y(t),`$`(t,2)) = -...

> initialdata(SimpSys);

TABLE([Pivots = [_C6 <> 0, _C2 <> 0, _C3-1 <> 0, _C...
TABLE([Pivots = [_C6 <> 0, _C2 <> 0, _C3-1 <> 0, _C...
TABLE([Pivots = [_C6 <> 0, _C2 <> 0, _C3-1 <> 0, _C...
TABLE([Pivots = [_C6 <> 0, _C2 <> 0, _C3-1 <> 0, _C...
TABLE([Pivots = [_C6 <> 0, _C2 <> 0, _C3-1 <> 0, _C...
TABLE([Pivots = [_C6 <> 0, _C2 <> 0, _C3-1 <> 0, _C...
TABLE([Pivots = [_C6 <> 0, _C2 <> 0, _C3-1 <> 0, _C...

Although there is no simplification as in the previous examples, rifsimp has found additional constraints that the initial conditions must satisfy ( Constraint and Pivots in the initialdata output), and obtained an equation for the time derivative of lambda(t) . Maple's dsolve[numeric] can then be used to obtain a numerical solution, in the manner of the previous example.

4. Lie Symmetry Determination (you need the newly updated ODEtools library for this example)

This example shows the use of rifsimp to assist in determination of Lie point symmetries of an ODE. Given the ODE:

> ODE:=2*x^2*y(x)*diff(y(x),x,x)-x^2*(diff(y(x),x)^2+1)+y(x)^2;

ODE := 2*x^2*y(x)*diff(y(x),`$`(x,2))-x^2*(diff(y(x...

We can obtain the system of PDE for it's Lie point symmetries via:

> sys := {coeffs(numer(DEtools[odepde](ODE, [xi,eta](x,y), y(x))),_y1)};

sys := {eta(x,y)*x*y^2+x^3*y*diff(eta(x,y),y)+eta(x...
sys := {eta(x,y)*x*y^2+x^3*y*diff(eta(x,y),y)+eta(x...
sys := {eta(x,y)*x*y^2+x^3*y*diff(eta(x,y),y)+eta(x...
sys := {eta(x,y)*x*y^2+x^3*y*diff(eta(x,y),y)+eta(x...
sys := {eta(x,y)*x*y^2+x^3*y*diff(eta(x,y),y)+eta(x...
sys := {eta(x,y)*x*y^2+x^3*y*diff(eta(x,y),y)+eta(x...
sys := {eta(x,y)*x*y^2+x^3*y*diff(eta(x,y),y)+eta(x...

We can then use rifsimp to greatly simplify this system:

> rifsys:=rifsimp(sys);

rifsys := TABLE([Solved = [diff(eta(x,y),`$`(x,2)) ...
rifsys := TABLE([Solved = [diff(eta(x,y),`$`(x,2)) ...
rifsys := TABLE([Solved = [diff(eta(x,y),`$`(x,2)) ...
rifsys := TABLE([Solved = [diff(eta(x,y),`$`(x,2)) ...

From this point, pdsolve can be used to obtain the infinitesimals for the Lie point symmetries via:

> pdsolve(rifsys['Solved']);

{xi(x,y) = -x*_C2+2*x*_C3-2*x*_C3*ln(x)+_C1*x+x*_C2...
{xi(x,y) = -x*_C2+2*x*_C3-2*x*_C3*ln(x)+_C1*x+x*_C2...

See Also

checkrank , initialdata , dsolve[numeric] , Rif , rifsimp[adv_options] , rifsimp[cases] , rifsimp[nonlinear] , rifsimp[options] , rifsimp[output] , rifsimp[overview] , rifsimp[ranking]