Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
ode solver help please

Subject: ode solver help please

From: Richard Crozier

Date: 7 Jan, 2010 11:58:20

Message: 1 of 3

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

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. 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?

Thanks,

Rich

Subject: ode solver help please

From: Steven Lord

Date: 7 Jan, 2010 14:40:27

Message: 2 of 3


"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 first-order 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; (V-iR)/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.soft-sys.matlab (CSSM) FAQ: http://matlabwiki.mathworks.com/MATLAB_FAQ

Subject: ode solver help please

From: Richard Crozier

Date: 7 Jan, 2010 15:14:25

Message: 3 of 3

Hi Steve, thanks for that very comprehensive reply! Unfortunately you state:

> > 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)
> >

The quantity psi is the flux linkage in a coil within the machine, the value of the flux linkage is dependent on the local magnetic field which is a function of various physical parameters of the machine, the position of the coil and the instantaneous value of the current ( i ) in the coil.

I can precompute values of psi for a range of position and current values using a specialist finite element package, but there is no explicit formula for this quantity and hence no derivative. The voltage (or emf in this case) is a direct result of the change in flux linkage and hence we cannot work backwards from the voltage to get psi.

If I were to provide a numerical gradient calculated from my lookup table for
d psi / dx, do you think this will seriosly affect the operation of the ode solver?

i.e. if I do something like

d psi / dt = ( d psi / d xR ) * ( d xR / dt )

              = gradient(psi values, suitable xR steps) * (vA - vT)

in place of an explicit formula, would this be acceptable or likely to cause serious problems?

p.s. I know I could do this in a simulink model, which I do have access to (but am not familiar with), but how does simulink deal with the same problem? i.e. how does it deal with such a situation behind the block diagram? Basically, is there an equivalent matlab function to the differentiator block in simulink?

Thanks again for your help!

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us