"Richard Crozier" <r.crozier@ed.ac.uk> wrote in message
news:hi4i8s$n46$1@fred.mathworks.com...
>I have a set of differential equations I wish to solve and I'm not exactly
>sure how to go about it and was wondering if anyone could offer some
>pointers, or good places to look for pointers or examples. Basically, I
>wish to simulate an electromechanical system comprising two parts coupled
>by magnetic forces. one of these parts is driven by a very stiff drive
>(i.e. I set the position and velocity of that part) and the other is
>attached to a surface via a spring. The two parts have a magnetic force
>linking them that varies with their relative position and velocity.
>
> The magnetic force is made up of two parts, Fm and Fe, one of which varies
> only with relative position (Fm), and the other which depends on both the
> position and the rate of change of another quantity which varies with
> position (Fe). I will refer to this quantity as psi(xR) where xR is the
> relative position of the two parts. Therefore we have:
>
> xA : position of sprung part relative to some absolute known position
> xT : position of driven part relative to the same absolute position
> vA = d xA / dt : velocity of sprung part vB = d xT / dt : velocity of
> driven part (applied by drive)
> xR = xT  xA : relative position
>
> Fm(xR) : magnetic force, can be drawn from a lookup table given xr
> psi(xR) : can also be drawn from a look up table given xr
>
> V = d psi / dt
You can rewrite this as:
d psi/dt = V
> This voltage is applied to an inductive circuit where R and L are known
> constants to find the current ( i ) on which the force depends:
>
> V  i R  L di / dt = 0
>
> so
> di / dt = ( V  i R ) / L
>
> Fe( i, xR ) : electrical force, a function of i and xR
>
> Fs = k xA : spring force, dependent on position and the spring constant,
> k
>
> The acceleration of the sprung part is then given by the sum of the forces
> divided by the mass (mA):
>
> d2xA / dt2 = (Fm + Fe + Fs) / mA
>
> The difficulty here, is that I have a quantity for which I do not have a
> derivative function, i.e I have psi at any position, but not a formula for
> d psi / dt.
Sure you do, as long as you know V or can compute it from other quantities
you know. And if you don't know V and can't compute it, then you don't have
a formula for di/dt that you can use to compute that quantity and I think
you're pretty well stuck.
> I can obtain
> d psi / dt from the chain rule so that:
>
> d psi / dt = ( d psi / d xR ) * ( d xR / dt ) = ( d psi / d xR ) *
> relative velocity
>
> = ( d psi / d xR ) * (vA  vT)
>
> but this still leaves the problem of getting d psi / d xR for one of the
> matlab ode solvers. I have some limited experience with the ode solvers
> such as ode45 etc. and only a little experience solving odes in general,
> and I am not sure what to do in this case. I have an inkling that a dde
> solver may possibly be what I'm after but I'm not sure. Can I perform my
> own (certainly inferior) numerical gradient calculation using past values
> inside my derivative function, or is this a very bad idea?
You have a system of three ODEs, correct?
d psi/dt = V
di / dt = ( V  i R ) / L
d^2xA / dt^2 = (Fm + Fe + Fs) / mA
Rewrite the last ODE as two firstorder ODEs:
d psi/dt = V
di / dt = ( V  i R ) / L
dxA/dt = (dxA/dt) % Trust me, this is useful
d (dxA/dt)/ dt = (Fm + Fe + Fs) / mA
Now let's define a vector Q with four elements:
Q = [psi; i; xA; dxA/dt]
Now your system becomes:
dQ/dt = [d psi/dt; di/dt; dxA/dt; d (dxA/dt)/dt]
dQ/dt = [V; (ViR)/L; Q(4); (Fm+Fe+Fs)/mA]
Q' = func(t, Q, other_quantities)
This looks a lot like the form that HELP ODE45 indicates it expects, doesn't
it? Now implement your function:
function Qprime = func(t, Q, ...
and use it in conjunction with ODE45. If you need further help, the
documentation for ODE45 (doc ode45) includes a number of examples that you
can use as a template.

Steve Lord
slord@mathworks.com
comp.softsys.matlab (CSSM) FAQ: http://matlabwiki.mathworks.com/MATLAB_FAQ
