Is there a numerical ODE solver facility in Maxima
Subject: Is there a numerical ODE solver facility in Maxima
From: Leo Butler
Date: Sat, 4 Apr 2009 12:05:48 +0100 (BST)
On Fri, 3 Apr 2009, Zach wrote:
< Ok, thanks everyone, this has been useful. I am still unsuccessful, but
< have a direction to work in. Jaime, thanks for the text, I don't know
< Portuguese but the examples still help. I talked with a friend who had some
< experience with Runge-Kutta stuff and he helped me out a bit.
<
< Here is the equation:
<
< diff( f(x), x, 2) = -k * sin( f(x) ) * cos( f(x) ) /* k is just some
< positive constant */
Zach, this equation is hamiltonian and can be re-written as
y'' = -k*sin(y)
where y=2*f.
(%i1) de1 : diff( f(x), x, 2) = -k * sin( f(x) ) * cos( f(x) );
2
d
(%o1) --- (f(x)) = - k cos(f(x)) sin(f(x))
2
dx
(%i2) de2 : subst(f(x)=y/2,de1);
2
d y y y
(%o2) --- (-) = - k cos(-) sin(-)
2 2 2 2
dx
(%i3) trigreduce(de2);
2
d y k sin(y)
(%o3) --- (-) = - --------
2 2 2
dx
In this form, it is the nonlinear pendulum
equation, which preserves energy
h = (y')^2/2 - k*cos(y).
(%i1) depends(y,t)$
(%i2) de : diff(y,t,2) = -k*sin(y)$
(%i3) v : diff(y,t)$
(%i4) h : v^2/2 - k*cos(y)$
(%i5) dh_dt : diff(h,t);
2
dy d y dy
(%o5) -- --- + k sin(y) --
dt 2 dt
dt
(%i6) subst(de,dh_dt);
(%o6) 0
< f'(-0.5) + K*sin(f(-0.5))*cos(f(-0.5)) = 0; /* K is just some other constant
< */
< -f'(0.5) + K*sin(f(0.5))*cos(f(0.5)) = 0;
Given the conservation of energy, you can view your BVP as a nonlinear
equation in the variables y_{+-}, v_{+-} (the angle and angular velocity
at the time +-0.5):
v_- + K*sin(y_-) = 0
-v_+ + K*sin(y_+) = 0
(v_+)^2 - k*cos(y_+) = (v_-)^2 - k*cos(y_-)
I would think that these equations will be easier to solve than the
original des.
-----
Note that you can solve your des in terms of elliptic functions, using
conservation of energy. Unfortunately, maxima is quite weak here, and
does not recognise the elliptic_f function that appears.
-----
Leo
<
< This leads me to believe that I can still use a shooting method, (provide an
< initial guess of f(-.5) which gives me f'(-.5), then play around with my
< guess of f(-0.5) until I get a solution that satisfies the BC at f(0.5)).
<
< Of course, a small angle approximation (sin(f(x)) * cos(f(x)) -> f(x)) can
< be made which makes this easy to solve (f(x) ~ cos(x)) and this can be used
< for my initial conditions for the shooting method.
<
< So, I think I know a way to proceed except I am not sure how to do it with
< rk. Maybe it can't be done? I guess I could always do exactly what I
< describe and just use an Euler integrator and take small steps in x... maybe
< the error won't be that bad.
<
< (I hope I made no mistakes above, I haven't gotten a result yet, so I am
< unsure that my thinking is 100% correct)
<
< Zach
<
< 2009/4/2 Martin Sch?necker <ms_usenet at gmx.de>
<
< > (%i1) ? rk;
< > -- Function: rk (<ODE>, <var>, <initial>, <domain>)
< > -- Function: rk ([<ODE1>,...,<ODEm>], [<v1>,...,<vm>],
< > [<init1>,...,<initm>], <domain>)
< > The first form solves numerically one first-order ordinary
< > differential equation, and the second form solves a system of m of
< > those equations, using the 4th order Runge-Kutta method. <var>
< > represents the dependent variable. <ODE> must be an expression
< > that depends only on the independent and dependent variables and
< > defines the derivative of the dependent variable with respect to
< > the independent variable.
< > The independent variable is specified with `domain', which must be
< > a list of four elements as, for instance:
< > [t, 0, 10, 0.1]
< > the first element of the list identifies the independent variable,
< > the second and third elements are the initial and final values for
< > that variable, and the last element sets the increments that
< > should be used within that interval.
< > If <m> equations are going to be solved, there should be <m>
< > dependent variables <v1>, <v2>, ..., <vm>. The initial values for
< > those variables will be <init1>, <init2>, ..., <initm>. There
< > will still be just one independent variable defined by `domain',
< > as in the previous case. <ODE1>, ..., <ODEm> are the expressions
< > that define the derivatives of each dependent variable in terms of
< > the independent variable. The only variables that may appear in
< > those expressions are the independent variable and any of the
< > dependent variables. It is important to give the derivatives
< > <ODE1>, ..., <ODEm> in the list in exactly the same order used for
< > the dependent variables; for instance, the third element in the
< > list will be interpreted as the derivative of the third dependent
< > variable.
< > The program will try to integrate the equations from the initial
< > value of the independent variable until its last value, using
< > constant increments. If at some step one of the dependent
< > variables takes an absolute value too large, the integration will
< > be interrupted at that point. The result will be a list with as
< > many elements as the number of iterations made. Each element in
< > the results list is itself another list with <m>+1 elements: the
< > value of the independent variable, followed by the values of the
< > dependent variables corresponding to that point.
< > There are also some inexact matches for `rk'.
< > Try `?? rk' to see them.
< > (%o1) true
< >
< > Best regards,
< > Martin
< >
< > _______________________________________________
< > Maxima mailing list
< > Maxima at math.utexas.edu
< > http://www.math.utexas.edu/mailman/listinfo/maxima
< >
<
--
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.