MATLAB Newsgroup

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

"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

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!

You can think of your watch list as threads that you have bookmarked.

You can add tags, authors, threads, and even search results to your watch list. This way you can easily keep track of topics that you're interested in. To view your watch list, click on the "My Newsreader" link.

To add items to your watch list, click the "add to watch list" link at the bottom of any page.

To add search criteria to your watch list, search for the desired term in the search box. Click on the "Add this search to my watch list" link on the search results page.

You can also add a tag to your watch list by searching for the tag with the directive "tag:tag_name" where tag_name is the name of the tag you would like to watch.

To add an author to your watch list, go to the author's profile page and click on the "Add this author to my watch list" link at the top of the page. You can also add an author to your watch list by going to a thread that the author has posted to and clicking on the "Add this author to my watch list" link. You will be notified whenever the author makes a post.

To add a thread to your watch list, go to the thread page and click the "Add this thread to my watch list" link at the top of the page.

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.

The newsgroups are a worldwide forum that is open to everyone. Newsgroups are used to discuss a huge range of topics, make announcements, and trade files.

Discussions are threaded, or grouped in a way that allows you to read a posted message and all of its replies in chronological order. This makes it easy to follow the thread of the conversation, and to see what’s already been said before you post your own reply or make a new posting.

Newsgroup content is distributed by servers hosted by various organizations on the Internet. Messages are exchanged and managed using open-standard protocols. No single entity “owns” the newsgroups.

There are thousands of newsgroups, each addressing a single topic or area of interest. The MATLAB Central Newsreader posts and displays messages in the comp.soft-sys.matlab newsgroup.

**MATLAB Central**

You can use the integrated newsreader at the MATLAB Central website to read and post messages in this newsgroup. MATLAB Central is hosted by MathWorks.

Messages posted through the MATLAB Central Newsreader are seen by everyone using the newsgroups, regardless of how they access the newsgroups. There are several advantages to using MATLAB Central.

**One Account**

Your MATLAB Central account is tied to your MathWorks Account for easy access.

**Use the Email Address of Your Choice**

The MATLAB Central Newsreader allows you to define an alternative email address as your posting address, avoiding clutter in your primary mailbox and reducing spam.

**Spam Control**

Most newsgroup spam is filtered out by the MATLAB Central Newsreader.

**Tagging**

Messages can be tagged with a relevant label by any signed-in user. Tags can be used as keywords to find particular files of interest, or as a way to categorize your bookmarked postings. You may choose to allow others to view your tags, and you can view or search others’ tags as well as those of the community at large. Tagging provides a way to see both the big trends and the smaller, more obscure ideas and applications.

**Watch lists**

Setting up watch lists allows you to be notified of updates made to postings selected by author, thread, or any search variable. Your watch list notifications can be sent by email (daily digest or immediate), displayed in My Newsreader, or sent via RSS feed.

- Use a newsreader through your school, employer, or internet service provider
- Pay for newsgroup access from a commercial provider
- Use Google Groups
- Mathforum.org provides a newsreader with access to the comp.soft sys.matlab newsgroup
- Run your own server. For typical instructions, see: http://www.slyck.com/ng.php?page=2