Solving ODE that uses a max() function

[ijustlovemath]

I have an ODE that I'd like to solve in labview, but it involves the max() function for the third DE. How can I incorporate this into the RK4 solver? I also will be applying this over a ten minute period (t here is in minutes), so would I put the RK4 solver inside a while loop with a shift register, or just set the end time to be 10+ the start time? I've included the table of variables + constants for your reference.

[johnsold]

With the max function in there it is no longer an "ordinary" differential equation.

Depending on how well behaved the solutions are, you might be able to use two equations - with and without the I2 term. Each time through the loop check the value of I2 and choose the appropriate equation for the next loop. If I2 tends to cross zero frequently this approach may not converge to a stable solution.

Lynn

[ijustlovemath]

They do tend to converge. So would I need a for loop to cover each time step? What would this look like in LV?

[johnsold]

What might work will probalby require some experimentation. I would start simple. Put the I2>0 test in the while loop you have now and see how things behave.  If that is not good enough then you may need to write your own Solver so you can look at each time step. The ODE Solver block diagram is available so you can see how they did things (some of the details are messy).

Lynn

[ijustlovemath]

Okay, so I've got progress. Apparently, the right side of my equation is formatted badly, but the Gmath library is so mangled I can't make heads or tails of it. Maybe a cursory look at my code would help nail down where I've gone wrong. I've attached it to this post.

[johnsold]

The ode ctl.ctl (typedef) file is missing.

What values are you using for that cluster? Without data it is hard to see what is going on.

Do you know what a valid solution is for those values?  What do the elements in output array represent?

Having all the panels and block diagrams set to maximize is rather annoying. When debugging a program like yours, I may have most of the panels and diagrams open simultaneously.  With all maximized it is hard to follow the flow.

Lynn

[ijustlovemath]

Sorry for missing that. The cluster contains four numerics, measurement, clock speed, A pump, B pump. The order should reflect that list. The elements in the output array should be a solution vector of (g,x,i1,i2). All I know about the solutions is that, for large clock speeds (it's the one wired to the eq comparison on the conditional terminal coming from the outside of the whole loop), g should tend to 100, and i1 should tend to 10. The other two have no physical interpretation so I'm unsure of how they work. Also, vi freezes (I think while running f(xt) picker) when step size is sufficiently low. Sorry for maximizing, I don't usually change window behavior from default.

[johnsold]

The f(x,t) picker.vi has at least two errors in the formulas. -1.3E-1(i1-10)+ is missing * as is the last formula.

However, fixing those does not change things.

What values do you use for Clock Speed, A pump and B pump. The pump values become part of the formula.

Your VI may never stop because testing floating point values for equality is problematic due to the way numbers are represented in binary. Numbers which you think should be equal may differ by a few parts in 10^-15 due to roundoff errors. The Clock Speed divided by 60 is likely to produce that sort of error if Clock Speed is not an integer multiple of 15. Dividing by 3 or 5 will result in an infinitely continuing fraction in binary.  One way to get around the issue is to subtact the two numbers, take the absolute value and then test for less than some small tolerance like 1E-15 or Machine Epsilon, found on the Numeric palette.

Lynn

[ijustlovemath]

Didn't think of that one! The floating point suggestion is valid. Clock speed is 600, a is 10 and b is 50. I set measurement to be 60. Can I just set the representation to double to fix the range error?

[johnsold · Accepted solution]

I do not understand your question about range error. Representations on all the numerics are already double.

I found the fundamental problem. Column to 1D array.vi returns an array with one element. That causes the RK solver to choke because it needs 4 elements. The fix is to just use Index Array. This version converges in about 11 iterations.

The solution seems to be sensitive to the value of h, so something is not quite right yet. Fixing that may be dependent on your knowledge of what the equations represent.

I am attaching the things I changed. In the f(x,t) picker.vi I added the * as mentioned previously. I also reformatted some of the constants, although I think that was unnecessary. In the Main VI I added a bunch of indicators so I could see what was happening. Those indicators led me to the array size I also added a Wait to allow me to see each iteration's results. You will want to remove that in your final version.  I implemented an "approximately equal" function to stop the loop. I also set it to stop on an error.

I would consider calculating the two possible formulas outside the loop and then just select the one you want inside. It does not take long to execute but if you needed hundreds or thousands of iterations to solve the equations it could add up.

Lynn