Differential-Algebraic equation help

Hello all,

I am brand new to using Igor Pro, so I know almost nothing about its interface and language. A paper that I read through disclosed a method for determining a particular parameter by running a simulation of a polymerization and fitting the simulated data to the experimental data by adjusting that parameter. The simulation involves solving a system of 6 ODEs with the Backwards Differentiation method, in addition to updating a number of algebraic equations at each time step during numerical integration. Could someone give me a brief tutorial on how I might go about solving the ODEs and algebraic equations simultaneously in Igor?

You have a long road ahead of you...

Start by reading the description of how to use IntegrateODE. Copy this command, paste it into Igor's command line and hit Enter:
DisplayHelpTopic "Solving Differential Equations"

If you haven't done it already, I strongly recommend selecting Help->Getting Started and doing at least the first half of the Guided Tour.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com
johnweeks wrote:
You have a long road ahead of you...

Start by reading the description of how to use IntegrateODE. Copy this command, paste it into Igor's command line and hit Enter:
DisplayHelpTopic "Solving Differential Equations"

If you haven't done it already, I strongly recommend selecting Help->Getting Started and doing at least the first half of the Guided Tour.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com


I already had read through the "Solving Differential Equations" tutorial and help section, but it never mentioned how to incorporate algebraic equations at each time step as well. Any quick tips?
The Derivative function embodies what you need to provide. The solver takes care of the time-stepping itself. By controlling the time steps, it can optimize the solution.

If you really have a system with derivatives that depend on the time of a given step, it is passed into your function via the second input parameter of the derivative function. In the examples, that is generally called either "xx" or "tt".

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com
johnweeks wrote:
The Derivative function embodies what you need to provide. The solver takes care of the time-stepping itself. By controlling the time steps, it can optimize the solution.

If you really have a system with derivatives that depend on the time of a given step, it is passed into your function via the second input parameter of the derivative function. In the examples, that is generally called either "xx" or "tt".

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com


I think I now know how to use the Derivative function to input my differential equations, but what I still am unsure of is how to incorporate non-differential equations into the same system and have those solved with the differential equations (the algebraic equations). Basically, some of the differential equations have terms in them that are dependent on other non-differential equations, so both the differential and algebraic equations needs to be solved simultaneously. For a short example, I might have:


dx/dt = Ax + y
dy/dt = By - z
z = C*x

So the first two are differentials, but the third is algebraic. I know in this example I could just substitute in the expression for z, but in my actual application the algebraic expressions are much more complex, so simplification is difficult. Therefore, I was wondering if Igor takes non=differential equations with the differentials, (i.e. could I list all three of those relationships in the Derivative function?)
My suspicion is that they can be implemented simply as equations assigned to a variable, then use the variable inside the expression for the derivative. It's hard to know without seeing the actual equations you are trying to solve.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com
johnweeks wrote:
My suspicion is that they can be implemented simply as equations assigned to a variable, then use the variable inside the expression for the derivative. It's hard to know without seeing the actual equations you are trying to solve.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com



I have attached here the equations in case that would help. I appreciate your time by the way; it would be so helpful if I could figure this thing out. Equations 9-14 are the differential equations, with the one non-integrable parameter being Tx and Ty. Tx and Ty can be computed from a series of algebraic relationships that could be simplified, which are equations 16, 18-19, 21-24. One more thing: although the differential equations need to be solved, it is actually the algebraic expression given in equation 24 that will yield the points I need for the plot I want to generate (the Disp variable vs. time). This equation is dependent on the differential equations, however, which is why I need to solve all of these in the first place.
equations.pdf (153 KB)
OK. So if I have understood correctly, you need to solve the ODEs in order to get the inputs to the algebraic equations. That will require making some additional waves and filling them via wave assignment expressions that include the relevant parts of the ODE solution. DisplayHelpTopic "Waveform Arithmetic and Assignment".

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com
johnweeks wrote:
OK. So if I have understood correctly, you need to solve the ODEs in order to get the inputs to the algebraic equations. That will require making some additional waves and filling them via wave assignment expressions that include the relevant parts of the ODE solution. DisplayHelpTopic "Waveform Arithmetic and Assignment".

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com


Thanks so much I'll give it a look.
johnweeks wrote:
OK. So if I have understood correctly, you need to solve the ODEs in order to get the inputs to the algebraic equations. That will require making some additional waves and filling them via wave assignment expressions that include the relevant parts of the ODE solution. DisplayHelpTopic "Waveform Arithmetic and Assignment".

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com


Hello again,

So I read through the waveform arithmetic and assignment section and I still have a few questions. Firstly, yes, the ODE solutions are inputs to algebraic equations, but the outputs of some of the algebraic equations are also inputs in the ODEs. For example, Tx and Ty are variables that are both in the ODEs as well as dependent on the ODEs through algebraic relationships. So I can't just solve the ODEs first, then plug in the outputs as inputs to the algebraic equations, because in order to solve the ODEs, I need to update the algebraic values along the way. So where should I algebraically define the waves Tx and Ty: within the derivative function or outside of it, but before it? Also, how can I perform operations on waves, such as multiplying two waves together? I can't find the proper syntax for this type of operation.
As an update to try to clarify, I posted a screenshot of the current derivative function I have created.

The highlighted "1" in the equation for dydt[2] is just a placeholder that is actually supposed to be a variable called Tx. This variable Tx is dependent on the differential state variables yw[1] and yw[2]. So the differential equation for dydt[2] requires the current value of Tx as an input, but Tx then requires the updated values of yw[1] and yw[2] as inputs for the next step. I am not sure how to incorporate this.
The inputs to your derivative function are a set of fixed parameters (in the examples, "pw"), the current value of the variable of integration ("xx" or "tt" in the examples) and the current value of the solution ("yw" in the examples). The output that you must compute is the set of derivatives ("dydx" or "dydt" in the examples). So I think the inputs to the derivative function have all the information you need.
Quote:
Also, how can I perform operations on waves, such as multiplying two waves together? I can't find the proper syntax for this type of operation.

If you have three waves, a, b, and c, you can fill c with the product of a and b with a simple assignment: c = a*b. That assumes that a, b, and c all have the same number of points. It also doesn't handle the case of multidimensional waves.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com
The input yw contains the current value of the solution. Your function will be called at various points in the state variable space, and the calls may jump back and forth in values of tt. I'm not sure what you mean by "Tx then requires the updated values of yw[1] and yw[2] as inputs for the next step". You get the current solution in yw; you can't have the solution to the next step until you get there.

It may be that you need a different method than we offer.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com
johnweeks wrote:
The input yw contains the current value of the solution. Your function will be called at various points in the state variable space, and the calls may jump back and forth in values of tt. I'm not sure what you mean by "Tx then requires the updated values of yw[1] and yw[2] as inputs for the next step". You get the current solution in yw; you can't have the solution to the next step until you get there.

It may be that you need a different method than we offer.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com


What i mean is this: the set of equations dydt[0-5] define the derivatives of the six inputs yw[0-5]. This is a system of ODEs with 6 ODEs. However, one of the terms in the derivative dydt[2] is not an input that has a differential equation like the other inputs yw[0-5], but rather an input that has an algebraic dependence on those other inputs. As a simplified version of the scenario:

dydt[0]=yw[0]*pw[0]
dydt[1]=yw[1]-yw[0]+Tx
Tx=yw[1]*yw[0]/pw[1]

So its a system of 2 ODEs, dydt[0] and dydt[1], but a third non-differential equation for Tx is necessary that has yw[0] and yw[1] as inputs.
Quote:

dydt[0]=yw[0]*pw[0]
dydt[1]=yw[1]-yw[0]+Tx
Tx=yw[1]*yw[0]/pw[1]

How about:
    Variable Tx=yw[1]*yw[0]/pw[1]
    dydt[0]=yw[0]*pw[0]
    dydt[1]=yw[1]-yw[0]+Tx


John Weeks
WaveMetrics, Inc.
support@wavemetrics.com
johnweeks wrote:
Quote:

dydt[0]=yw[0]*pw[0]
dydt[1]=yw[1]-yw[0]+Tx
Tx=yw[1]*yw[0]/pw[1]

How about:
    Variable Tx=yw[1]*yw[0]/pw[1]
    dydt[0]=yw[0]*pw[0]
    dydt[1]=yw[1]-yw[0]+Tx


John Weeks
WaveMetrics, Inc.
support@wavemetrics.com


Perfect thanks. And I can define this variable within the derivative function? Lastly, will the changing values of Tx be stored in an array so that I can view the entire range of Tx values over time?
integratemate wrote:
And I can define this variable within the derivative function?

Yes. Please read the Programming help file.
Quote:
Lastly, will the changing values of Tx be stored in an array so that I can view the entire range of Tx values over time?

No. Only the solution values will be saved. But since you have an expression for generating Tx from the solution values, you can re-generate those values after IntegrateODE finishes.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com
johnweeks wrote:
integratemate wrote:
And I can define this variable within the derivative function?

Yes. Please read the Programming help file.
Quote:
Lastly, will the changing values of Tx be stored in an array so that I can view the entire range of Tx values over time?

No. Only the solution values will be saved. But since you have an expression for generating Tx from the solution values, you can re-generate those values after IntegrateODE finishes.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com



Alright the procedure with all of my new variables compiles fine, but now the IntegrateODE command results in an error saying the step size has become too small. I tried raising the error tolerance significantly with no success. Any final tips? I know this exact set of equations has been done in Igor before because it was in the experimental section of another paper, but the detail is very limited so I don't know how they did it. Here is what I have in the procedure that compiles fine. Is there any glaring error that you see that would make you think this can't be done (see attached)