Alvaro, and interested others,

Here are the promised files to see if there are SMATH alternates to the coupled finite difference equations used in each example.

First is a calculation of the yield of a fission weapon. To put everyone's mind at ease, there is nothing classified. Everything is available in open literature. My model was derived from the math in The Physics of the Manhattan Project by Bruce Cameron Reed. He was very kind to correspond with me when I had specific questions. I wrote this specific sm file to estimate what yield Iran might be able to produce if they develop compression similar to that of the TRINITY Test. The results are disconcerting.

This file runs in about 27 seconds on my 10 year-old PC.

FISSION YIELD WITH COMPRESSION.sm (148 KiB) downloaded 70 time(s).

This file calculates the High Altitude EMP (HEMP) coupling to a cable. Again, all is available in the open literature. NOTE: The HEMP fields need correction based on more recent data. I will include that when I get the time. For now, this serves as an example of the time domain finite difference solution to the coupled transmission line equations. This is the most difficult example because the coupling is in both time and position on the line.

This file runs in about 19 seconds on my PC.

HEMP TDFD CABLE COUPLING.sm (167 KiB) downloaded 71 time(s).

This final example was inspired during the COVID-19 pandemic. It is an SIR epidemiology model, where S=susceptible population, I=infected population, and R=the removed population by either recovery or death. The recovery time is 14 days. I removed the extrapolation to future dates for simplicity. Also, confidence was low projecting into the future because the virus reproduction rate, Ro, was constantly changing.

This file runs in about 2-3 seconds on my PC.

SIR EPIDEMIOLOGY FOR COUPLED DE.sm (923 KiB) downloaded 58 time(s).

For anyone interested, I have created mov presentations of these three issues plus my climate studies. Here is a link to my folder on Microsoft's OneDrive.

https://1drv.ms/f/c/86bd4c6ece5074bb/Ert0UM5uTL0ggIZrAAAAAAABSZS0RIW3H_jmVhSU_tPfjQ

The studies make use of my MATHCAD-15 files from which I developed these SMATH files. I am gradually converting to SMATH because PTC refuses to help me get MATHCAD-15 on a new computer if this one fails, even though I have a perpetual license. At present, there is more detail in my MATHCAD files than the SMATH files, but I am getting there.

Regard,
Reg

Alvaro,
Here are the equations for the example problems.

Equations for Example DE Problems .pdf (1 MiB) downloaded 76 time(s).

Let me know if you need more.

Reg

Hi Reg. Thanks! This would be the first approximation to the covid ode. It is a family of 57 coupled differential equations, there are surely several ways to approach the problem, but at least here it is shown that it is possible to solve coupled odes with SMath. For the fission problem I need more time, and for the transmission lines problem I probably need more lives.



SIR EPIDEMIOLOGY FOR COUPLED DE - RK.pdf (432 KiB) downloaded 65 time(s).

SIR EPIDEMIOLOGY FOR COUPLED DE - RK.sm (1 MiB) downloaded 64 time(s).

Best regards.
Alvaro.

Hi. I just extract the parameters, and try to solve the ode. But somethings goes wrong, and can't see why.

FISSION YIELD WITH COMPRESSION just the ode.sm (167 KiB) downloaded 54 time(s).



Best regards.
Alvaro.

Wrote

... I do have a question. On the SIR model, is the coupled ODE approach faster than the finite difference model? ...

It is quite a bit slower. The "dsol" method involves solving a system of non-linear equations with nleqsolve and then calling the ode solver. However, in this case there would be no need to use it, since the derivatives of the functions are already "cleared", so the ode solver can be used directly, and in this case it is much faster.



SIR EPIDEMIOLOGY FOR COUPLED DE - RKFIXED.sm (1 MiB) downloaded 75 time(s).

Best regards.
Alvaro.

Alvaro,

Thanks much.

Two final questions.

1. How do you decide between rkfixed and rkadapt?

2. I don't understand this construction for beta (t).



Is there an alternate way to define to define beta(t)?

Thanks again.
Reg

I am still trying to understand why Alvaro's ODE solver for the SIR coupled ODEs does not give the same result as my finite difference method (contained in the same sm file above). I must admit that I almost never use ODE solvers (in MATHCAD or SMATH); therefore, I don't understand them as well as finite difference.

So, I decided to start with a simpler problem where the answer is known. This sm file was derived based on the book, Smath for Physics, A primer. The problem is referenced inside the sm file. The book is also referenced on the SMath Studio main webpage.

COUPLED ODEs FROM SMATH FOR PHYSICS.sm (40 KiB) downloaded 63 time(s).

This gives the correct answer. So, I continue my quandary of why Alvaro's coupled ODE model does not fit the data like the finite difference method. As far as I can tell with my infantile understanding of SMATH ODE solvers, they should give the same result.

If anyone knows the answer, I would appreciate the answer.

Finally, does SMATH run well on Windows 11?

Regards,
Reg

Wrote

...1. How do you decide between rkfixed and rkadapt?

2. I don't understand this construction for beta (t). Is there an alternate way to define to define beta(t)? ...

Hi Reg.

1. The author of these functions, Viacheslav, once commented that he wanted to reproduce the same values ​​that would be obtained with the Mathcad functions. Besides those two, he has about 5 more plugins with other functions to solve differential equations. In all of them the functions always return N intervals for the independent variable, that is, N+1 points equally distributed, so if there was any interesting information about the original problem in the independent variable, regarding where more points were used and where less, this is lost, and it does not serve as a criterion to choose between adaptive or fixed-step methods. One detail about rkadapt is that sometimes it fails or "gets stuck", while rkfixed is more "stable", although many times the results with rkadapt are more precise. Hence my highly scientific and specialized criterion is the following: if there are no visually appreciable differences between the graph with rkfixed and rkadapt, I choose rkfixed. And if it is very slow, I try with "Adams", which is usually much faster. As for stability, convergence, and other studies, I don't usually do them. And I suspect that maybe the solution is not the right one, I resort to some function in the other Viacheslav plugins.

2. There are other ways, but almost all of them involve the use of "if", and those ways cannot be used, since SMath's native "if" cannot be passed to the other functions. Other ways could use findrows, or max, vectorizing t
Best regards.
Alvaro.

Hi Alvaro,

Thanks for the additional information.

I posted a simpler comparison of finite difference vs the ODE solver for the SIR problem. It’s posted in the Questions section.

To keep it simple, I made Ro (beta) constant and I omitted any data. It just compares the two methods. I don’t understand why they don’t agree because the ODE solver exactly solves the problem from the book SMath for Physics that I posted above in this thread.

I would really like to switch to the ODE solver for most problems because it is much faster…about a factor of ten in the comparison in the Questions section. It’s very obvious if you really increase the number of steps to around 5000. The finite difference takes more than 40 seconds while the RK runs in less than 4 seconds. I hope there is something simple to fix the inconsistency.

Thanks again.
Reg

Hi Reg. It turns out that one can use "if" with rkfixed. It seems that the trick is to use "line" in the definition of the D function. This version should also work in mathcad. In addition, a very simple routine of the RK23 method is added.



SIR EPIDEMIOLOGY FOR COUPLED DE - RK23.sm (1 MiB) downloaded 64 time(s).

Best regards.
Alvaro.

Wrote

I am still trying to understand why Alvaro's ODE solver for the SIR coupled ODEs does not give the same result as my finite difference method ...

Hi Reg. My guess is that the Ro values in SEGMENTS matrix would be specific to the finite difference method and would not be accurate if an RK method is used.

Best regards.
Alvaro.

Hi Alvaro,
Thanks. I agree.

I am encouraged that it works well for the Fission and Climate problems.

I think I now have enough to keep me busy until I exit the planet. I have a new Windows 11 PC on order. I dread all the work to get all my apps transferred over and getting it configured the way I want it. It seems like when I finally get a computer configured the way I want it, it is close to becoming obsolete.

Take care and thanks again.
Reg

Wrote

Wrote

I am still trying to understand why Alvaro's ODE solver for the SIR coupled ODEs does not give the same result as my finite difference method ...

Hi Reg. My guess is that the Ro values in SEGMENTS matrix would be specific to the finite difference method and would not be accurate if an RK method is used.

Best regards.
Alvaro.

I guess that the update frequency of the rate equations is part of the model just like in financial calculations with interest rates. The frequency of accumulation of a given nominal interest affects the effective interest rate and probably the model parameters can be converted into equivalent values for different update frequencies. In numerics you usually try to get the time step converged, so that the results aren't affected by this choice.

Hi Martin,

This is a problem of fitting a model to data and not just a math problem. Since the data was only available on a daily basis and later CDC only gave it weekly and then stopped altogether, the step is one day at best.

However, the real answer to the difference is that the SEGMENTS data does not get into the ODE formalism correctly. Someone else like Alvaro will need to explain that because I tend to use FD for all coupled differential equations. If they are coupled in both time and space, as in the EMP example, I don't think there is an ODE formalism that exists to handle that. But, I am much less familiar with ODE solvers. I don't even use them in MATHCAD.

Nevertheless, my example below demonstrates that the SEGMENTS data is not getting into the ODE formalism. Look down through the file and you will see that I have arbitrarily changed the final Ro to an absurd value of 20. Now, even the TDFD does not fit the last data point. But there is no change to the ODE solution; therefore, it's not getting the change and thus it cannot fit the data. Ro in each segment was adjusted to get a good fit to the data over the data range provided, usually days or weeks. If you put the best fit value back into the last Ro, the fit goes back to an excellent fit at the last data point. Further down in the file, you can see the actual difference between the model and the data. With the best fit Ro, the difference is about 8E-5, clearly less than a person. If Ro = 20, the difference is about 37 million people.

In either case, the ODE model does not change. See the purple box at the bottom as you select the final Ro value. The TDFD method will change accordingly.

As I understand from Alvaro, it's very difficult to define a changing parameter outside the ODE block and have the ODE block correctly get that parameter; however, as I said, I seldom use ODE solvers. Someone else will need to explain that.

SIR COUPLED TDFD VS ODE.sm (937 KiB) downloaded 67 time(s).

Regards,
Reg