This is an outgrowth of the thread under Samples, al_nleqsolve issues. That thread evolved into a discussion of finite difference methods versus the SMATH ODE solvers. I have simplified the coupled differential equations so there is no nesting in the finite difference approach and the reproduction rate, Ro, is constant. This file has three coupled differential equations solved by finite difference and the ODE Runge-Kutta solver. They do not agree. Does anyone see the cause of the difference?

SIR FINITE DIFFERENCE vs ODE.sm (57 КиБ) скачан 35 раз(а).

Thanks,
Reg

Try using a smaller timestep. Image below shows comparison with a timestep of 1/10 day as opposed to your original 1 day.

Note that your "Finite Difference" method is the simple Euler method, which is generally considered to be inaccurate when the curves are changing rapidly. A smaller timestep often helps, but that can only be taken so far before numerical precision limits further improvement.

Note also that the RK4 method is also a finite difference method (a higher order one).




SIR FINITE DIFFERENCE vs ODE_b.sm (57 КиБ) скачан 33 раз(а).

Thanks, but why did you have to redefine my finite difference equations? If I just make my time step 1/10 and increase the steps to 1000, I don't get your result. So I am still confused.

Nevertheless, my ultimate goal is to get SMATH ODE to solve the following SIR, FUSION, and EMP problems.

Here is the SIR problem. It needs to have Ro changing throughout the calculation. Ro changes for specific ranges of dates to fit the actual CDC data. It only remains constant through each specific range.

SIR FINITE DIFFERENCE vs ODE.sm (56 КиБ) скачан 30 раз(а).

So far, there is no SMATH ODE solver that matches the data like the finite difference method. See the thread under Samples re Resonar.


Here is the FUSION problem.
FISSION YIELD WITH COMPRESSION.sm (148 КиБ) скачан 34 раз(а).

And here is the EMP problem.
HEMP TDFD CABLE COUPLING.sm (167 КиБ) скачан 26 раз(а).

I would be happy just to find an SMATH ODE solver that can match the data in the SIR problem like the finite difference method.

Anyway, I really appreciate your response. I would stick with MATHCAD except for one problem. I will shortly have to purchase a new Windows 11 computer and PTC will not help me get MATHCAD 15 installed on a new computer, even though I have a perpetual license. They claim there is some legal issue. They have a history of poor support. While the finite difference solutions for the above problems run very quickly on MATHCAD with significantly finer steps, they are extremely slow on SMATH. Thus my objective to find an ODE solver around this problem. In my above FUSION and EMP SMATH problem examples, I have sacrificed accuracy in to get a reasonable runtime in SMATH.

Thanks much.
Reg

Great! Thanks much.

Now, do you have a suggestion on doing the full SIR calculation with the ODE solver including the SEGMENTS data for Ro.

Thanks again. I appreciate you help.
Reg

Hi Reg. The fusion problem.



FISSION YIELD WITH COMPRESSION.sm (274 КиБ) скачан 27 раз(а).

Best regards.
Alvaro.

Alvaro,
Great. That is very close to my result if I increase the core expansion steps from 50 to 100. My runtime increases from about 27 seconds to about 45 seconds.



The MATHCAD result is here for 5000 steps in the core expansion.



Trying to use 5000 core expansion steps really would stress the SMATH runtime. That is the reason I keep looking for an ODE solver in SMATH. If I disable my finite difference calculation in your code it runs for about 20 seconds for 50 steps and 32 seconds for 100 steps, a significant improvement over the finite difference method in SMATH. Plus, it's probably more accurate.

Increasing the core expansion steps in the SMATH version reduces the accuracy of the linear fit to alpha and lambda; this does not happen in my MATHCAD version. If you look at the solve function that I use in SMATH for the fit, there is probably a better fit function in SMATH, but I was not very familiar with the fitting functions in SMATH when I developed this fission model. The complicated expression that obtains alpha comes from a solution of the neutron diffusion equation in spherical coordinates. In MATHCAD, I don't have to break up the f(alpha) transcendental equation into three components.

This is how it is done in MATHCAD.


Now, I have a race in time to determine if I reach the pine box after I understand you ODE method or not.

Great work and thanks as always.
Reg