coupled differential equation and Matlab ode45

Hi

I am solving 2 coupled ode's:

dKdt=P*(K0-K)/(1-t)

dG/dt=Q*(G0-K)/(1-t)

P, Q, K0, G0 are constants.

I am using the following functions:

function ISOdem(asp,KK)
variable asp
wave KK
Make/D/O/N=(10001,2) DEMout
SetScale/P x 0,0.0001,DEMout
SetDimLabel 1,0,K,DEMout       
SetDimLabel 1,1,G,DEMout                           
kk[4]=asp
DEMout[0][%K] = kk[0]          
DEMout[0][%G] = kk[1]                                                          
IntegrateODE/M=0/u=1000 Berryman, KK, DEMout
end

function Berryman(pw, tt, yw, dydt)
    Wave pw                    
    Variable tt                
    Wave yw                    
    Wave dydt                  
    variable asp=pw[4]
    variable theta, fn, nu, r, a, b, Pa, Qa

// lots of calculations to get P and Q

    theta=(asp/((1-asp^2)^(3/2)))*(acos(asp)-asp*sqrt(1-asp^2))
    fn=(asp^2/(1-asp^2))*(3*theta-2)
    nu=(3*pw[0]-2*pw[1])/(2*(3*pw[0]+pw[1]))
    r=(1-2*nu)/(2*(1-nu))
    a=pw[3]/pw[1]-1  // pw[3] = 0 so it is always true that a=-1
    b=(1/3)*(pw[2]/pw[0]-pw[3]/pw[1]) // pw[3] = 0 so it is always true that b=(1/3)*(pw[2]/pw[0])
       // also, when pw[2] = 0 (b = 0) there is virtually no discrepancy

    variable f1a=1+a*((3/2)*(fn+theta)-r*((3/2)*fn+(5/2)*theta-(4/3)))
    variable f2a=1+a*(1+(3/2)*(fn+theta)-(r/2)*(3*fn+5*theta))+b*(3-4*r)+(a/2)*(a+3*b)*(3-4*r)*(fn+theta-r*(fn-theta+2*theta^2))
    variable f3a=1+a*(1-(fn+(3/2)*theta)+r*(fn+theta))
    variable f4a=1+(a/4)*(fn+3*theta-r*(fn-theta))
    variable f5a=a*(-fn+r*(fn+theta-(4/3)))+b*theta*(3-4*r)
    variable f6a=1+a*(1+fn-r*(fn+theta))+b*(1-theta)*(3-4*r)
    variable f7a=2+(a/4)*(3*fn+9*theta-r*(3*fn+5*theta))+b*theta*(3-4*r)
    variable f8a=a*(1-2*r+(fn/2)*(r-1)+(theta/2)*(5*r-3))+b*(1-theta)*(3-4*r)
    variable f9a=a*((r-1)*fn-r*theta) + b*theta*(3-4*r)

    Pa = f1a/f2a
    Qa = 1/5*(2/f3a+1/f4a+(f4a*f5a+f6a*f7a-f8a*f9a)/f2a/f4a)

// the differential equation:
    dydt[0] = (pw[2]-yw[0])/(1-tt)*Pa
    dydt[1] = (pw[3]-yw[1])/(1-tt)*Qa
    return 0
end

 

It seems to be working fine, but when I compare the output with results from Matlab (its original environment) there is a discrepancy. This issue only seems to arise for certain cases (when 'asp' is small AND pw[2] is nonzero), which really should only affect the calculate values of P and Q, but otherwise not change how the function operates.

Is there any known difference between ODE solvers? Both Igor and Matlab are supposedly using Runge-Kutta (plus, I get the same result from Igor using other algorithms).

Thanks!

If you find that the Igor solution matches Matlab mostly, but has discrepancies with certain conditions, it would seem reasonable to look for errors in your code under those conditions.

Assuming there are no errors in the code, you might also want to look at the control of truncation error (DisplayHelpTopic "Error Monitoring"). By default, Igor sets the error tolerance to 1e-6, and constant error scaling equal to 1. When using /M=0 or /M=1 the error for a given step that is compared to the error limit is the maximum of the error estimates for each of the components. Using constant error scaling means that when the solution is small (close to zero) the errors as a percentage will be larger than when the solution is large.

The treatment of truncation error varies depending on the solver you use. The RK solver in Igor is based on Numerical Recipes in C, and it is from that book that we got the max(error_i) technique. I think (though I am not certain by any means) that Matlab uses the SUNDIALS library. That library use the RMS error. And by default, the truncation error combines scaling with a constant (absolute error) with scaling by the magnitude of the solution (relative error). That may keep a tighter tolerance when the solution gets small.

You can use the /F flag to add a relative tolerance term to IntegrateODE, and the /S flag can specify a wave to use for scaling the parts of the error.

If you have no coding errors, then I would look at the truncation error control as the source of your problem.

Thanks for the quick response.

I just found the source of the problem: P and Q are not supposed to be constant (as I set them), but are functions of the yw[0] and yw[1] I am solving for. In other words, the coupled equations are:

dydt[0] = (pw[2]-yw[0])/(1-tt)*Pa(yw[0])
dydt[1] = (pw[3]-yw[1])/(1-tt)*Qa(yw[1])

I am not sure of how to incorporate this.

Can I just turn all the variables (theta, fn, nu, r, a, b, Pa, Qa) into waves and replace pw[0] and pw[1] with yw[][0] and yw[][1]?

 

Thanks. 

I will formulate my question carefully tomorrow and include my new script if I can’t get it working - but I think my proposal is what you are suggesting (using the 2 columns of yw in the expressions for Pa and Qa - which originally were variables but actually need to be waves).

Thanks again

 

 

 

 

I'm glad it was an easy fix!

FYI, I wasn't understanding how yw[0] and yw[1] were used in the function, so I thought I had to make all the expressions involving yw into waves, instead of keeping the variables (of so-called "common subexpressions") as they are and just inserting yw[0] or yw[1] into the places they belong (i.e., in place of pw[0] or pw[1]).

Thanks for your help, and the help of the Solving Differential Equations your team has produced.

Yes, the input yw wave for the derivative function is confusing. During the integration, Igor actually creates a new wave with as many rows as your output wave has columns (read that carefully!). Every time it calls your derivative function, it sets the wave elements to whatever point in your solution space is being worked on. That isn't necessarily anything in your solution wave, and won't necessarily even wind up in your solution wave. But as the integration moves forward, it must compute results at a variety of points between where it was and where it is going. Those points are the values put into yw.

Especially with the methods other than RK, but also with RK to some extent, the values aren't necessarily even in order in terms of the X variable. You really can't assume anything about where the derivative function will go next.