coupled differential equation and Matlab ode45

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

The integrand function has yw as an input. Just use that.

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.