ode
-- the domain of ordinary
differential equationsode
(eq, y(x))
represents an ordinary
differential equation (ODE) for the function y(x)
.
ode
({eq1, eq2, ...}, {y1(x), y2(x), ...})
represents a system of ODEs for the functions y1(x)
,
y2(x)
etc.
ode(eq, y(x))
ode({eq <, inits>}, y(x))
ode({eq1, eq2, ... <, inits>}, {y1(x), y2(x),
...})
eq, eq1, eq2, ... |
- | equations or arithmetical
expressions in the unknown functions and their derivatives with
respect to x . An arithmetical expression is regarded as an
equation with vanishing right hand side. |
y, y1, y2, ... |
- | the unknown functions: identifiers |
x |
- | the independent variable: an identifier |
inits |
- | the initial or boundary conditions: a sequence of equations |
an object of type ode
.
numeric::odesolve
, numeric::odesolve2
, plot::ode
eq
, eq1
etc., the
unknown functions must be represented by y(x)
,
y1(x)
etc. Derivatives may be represented either by the
diff
function or by the
differential operator D
.
Note that the token '
provides a handy short cut:
y'(x) = D(y)(x)
= diff(y(x),
x)
.x
. Multivariate expressions such as y(x,
t)
are not accepted.D
(or the token
'
) must be used to specify values of derivatives at some
point. E.g.,
y(1) = 2, y'(0) = 0, y''(0) = 1
is a valid sequence of boundary conditions for inits
.
Boundary conditions of the first and second kind are allowed. Mixed conditions are not accepted.
The initial/boundary points and the corresponding initial/boundary values may be symbolic expressions.
ode(
{eq, inits}, y(x))
to specify the
conditions.ode
domain is to provide an
environment for overloading the function solve
.
In the case of one single equation (possibly together with initial
or boundary conditions), solve
returns a set of explicit
solutions or an implicit solution. Each element of the set represents a
solution branch.
In the case of a system of equations, solve
returns a set of lists of
equations for the unknown functions. Each list represents a solution
branch.
An symbolic solve
call is returned if no solution is found.
setuserinfo(ode, 10)
, a solve
command provides information on
MuPAD's way of solving ODEs.In the following, we show how to create and solve a
scalar ODE. First, we define the ODE x^2*y'(x) + 3*x*y(x) =
sin(x)/x. We use the quote token '
to represent
derivatives:
>> eq := ode(x^2*y'(x) + 3*x*y(x) = sin(x)/x, y(x))
/ sin(x) 2 \ ode| 3 x y(x) - ------ + x diff(y(x), x), y(x) | \ x /
We get an element of the domain ode
which
we can now solve:
>> solve(eq)
{ C1 + cos(x) } { - ----------- } { 3 } { x }
>> delete eq:
An initial value problem is defined as a set consisting of the ODE and the initial conditions:
>> ivp := ode({f''(t) + 4*f(t) = sin(2*t), f(0) = a, f'(0) = b}, f(t))
ode({f(0) = a, D(f)(0) = b, 4 f(t) + diff(f(t), t, t) - sin(2 t)}, f(t))
>> solve(ivp)
{ sin(2 t) b sin(2 t) t cos(2 t) } { -------- + a cos(2 t) + ---------- - ---------- } { 8 2 4 }
>> delete ivp:
With some restrictions, it is also possible to solve systems of ODEs. First, we define a system:
>> sys := {x'(t) - x(t) + y(t) = 0, y'(t) - x(t) - y(t) = 0}
{y(t) - x(t) + D(x)(t) = 0, D(y)(t) - y(t) - x(t) = 0}
A call to solve
yields the general solution with
arbitrary parameters:
>> solution := solve(ode(sys, {x(t), y(t)}))
{[y(t) = C7 exp((1 + I) t) + C8 exp((1 - I) t), x(t) = I C7 exp((1 + I) t) - I C8 exp((1 - I) t)]}
To verify the result, we substitute it back into the
system sys
. However, for the substitution, it is necessary
to rewrite the system into a notation using the diff
function:
>> eval(subs(rewrite(sys, diff), op(solution)))
{0 = 0}
>> delete sys, solution:
In this example, we point out the various return formats
of ode
's solve facility. First, we solve an ODE with an
initial condition. The solution involves a symbolic integral:
>> solve(ode({y'(x) + x*y(x) = cos(x), y(0) = 3}, y(x)))
{ / 2 \ 2 1/2 } { | x | int(cos(t2) exp(t2 ) , t2 = 0..x) } { 3 exp| - -- | + ----------------------------------- } { \ 2 / 2 1/2 } { exp(x ) }
The following system is solved incompletely:
>> sys := {x'(t) = -3*y(t)*z(t), y'(t) = 3*x(t)*z(t), z'(t) = -x(t)*y(t)}: solution := solve(ode(sys, {x(t), y(t), z(t)}))
2 1/2 2 1/2 {[x(t) = (C10 - C11 + 3 z(t) ) , y(t) = (- C10 - 3 z(t) ) ] 2 1/2 , [x(t) = (C10 - C11 + 3 z(t) ) , 2 1/2 y(t) = - (- C10 - 3 z(t) ) ], 2 1/2 [x(t) = - (C10 - C11 + 3 z(t) ) , 2 1/2 y(t) = (- C10 - 3 z(t) ) ], [ 2 1/2 x(t) = - (C10 - C11 + 3 z(t) ) , 2 1/2 y(t) = - (- C10 - 3 z(t) ) ]}
In these four different solutions branches, no solution
for the unknown function z(t)
is provided. In fact, the
partial solution above is subject to a further condition on
z(t)
. We substitute the first of our four solutions back
into the system sys
:
>> eval(subs(rewrite(sys, diff), solution[1]))
{ 2 1/2 { diff(z(t), t) = - (- C10 - 3 z(t) ) { { 2 1/2 3 z(t) diff(z(t), t) (C10 - C11 + 3 z(t) ) , - -------------------- = 2 1/2 (- C10 - 3 z(t) ) 2 1/2 3 z(t) (C10 - C11 + 3 z(t) ) , 3 z(t) diff(z(t), t) 2 1/2 } ------------------------ = - 3 z(t) (- C10 - 3 z(t) ) } 2 1/2 } (C10 - C11 + 3 z(t) ) }
For each solution branch, there remains is exactly one
differential equation for z(t)
to solve.
>> delete sys, solution:
It may happen that MuPAD cannot solve a given
equation. In such a case, a symbolic solve
command is returned:
>> solve(ode(y'(x) + y(x)^2 = b + a*x, y(x)))
2 solve(ode(- b - a x + diff(y(x), x) + y(x) , y(x)))
MuPAD's ODE solver contains algebraic algorithms for computing Liouvillian solutions of linear ordinary differential equations over the rational functions as well. These algorithms are based on differential Galois theory. For second order equations, an algorithm similar to the Kovacic algorithm is implemented. However, instead of computing the so-called Riccati-polynomial of a solution, this algorithm computes the solutions directly using formulas. Only when both solutions are algebraic functions, then in three cases it is necessary to compute the minimal polynomial of a solution, again using formulas. Hence, for Kovacic's famous example
y''(x) + (3/(16*x^2) + 2/(9*(x - 1)^2) - 3/(16*x*(x - 1)))*y(x) = 0,
it is possible to compute the minimal polynomial of solution:
>> solve(ode(y''(x) + (3/(16*x^2) + 2/(9*(x - 1)^2) - 3/(16*x*(x - 1)))*y(x), y(x)))
24 8 8 {C2 RootOf(_Y1 + 2985984 x (x - 1) - 4 8 6 2 16 3 2799360 x _Y1 (x - 1) - 4320 x _Y1 (x - 1) - 5 4 1/2 1/2 7 165888 x _Y1 (x - 2) (I x 3 - I 3 ) + 3 12 1/2 1/2 1/2 4 5760 I x _Y1 3 (x - 2) (I x 3 - I 3 ) , _Y1)}
MuPAD may find Liouvillian solutions for higher order equations as well. However, there is no guarantee that all of them are found. The following third order equation can be solved completely:
>> solve(ode(diff(y(x), x, x, x) + diff(y(x), x, x)*(2*x^2 - 1)/(2*x^2 - 2*x) + diff(y(x), x)*(295*x - 491*x^2 + 196*x^3 - 98) /(196*x^2 - 392*x^3 + 196*x^4) + y(x)*(3*x - x^2 - 1) /(196*x^2 - 392*x^3 + 196*x^4), y(x)))
{ { 1/2 1/14 { C3 (2 x + 2 (x (x - 1)) - 1) + { { { { C4 -------------------------------- + 1/2 1/14 (2 x + 2 (x (x - 1)) - 1) / | 1/2 1/14 | C2 (2 x + 2 (x (x - 1)) - 1) | \ / | 1/2 | int(exp(-x) (x (x - 1)) | \ 1/2 1/14 (2 x + 2 (x (x - 1)) - 1) , x) / 1/2 1/7 (2 x + 2 (x (x - 1)) - 1) - / 1/2 \ \ \ } | exp(-x) (x (x - 1)) | | | } int| --------------------------------, x | | | / 28 } | 1/2 1/14 | | | } \ (2 x + 2 (x (x - 1)) - 1) / / / } } }