Modeling uncertainty from known physics

Question

I would like to request some help from the data science community.

My task with machine learning is do achieve the following in pyTorch:

I have an equation given by:

$$ \frac{\mathrm{d} s}{\mathrm{d} t}=4a−2s+\lambda(s) $$

Where $a$ is an input constant and $\lambda$ is a non-linear term that depends on $s$.

I know that the true solution for $\lambda\left(s\right)$ is $\sin \left( s\right) \cos\left(s\right)$

I have generated data from [0, 5]s with the true solution and with the equation excluding the non-linear term for a range of initial conditions and a range of $a$'s.

import numpy as np
import pandas as pd
from scipy.integrate import solve_ivp

def foo_true(t, s, a):
    ds_dt = 4*a -2*s + np.sin(s)*np.cos(s)
    return ds_dt

def foo(t, s, a):
    ds_dt = 4*a -2*s
    return ds_dt


# Settings:
t = np.linspace(0, 5, 1000) 

a_s = np.arange(1, 10)
s0_s = np.arange(1, 10)

# Store the data
df_true = pd.DataFrame({'time': t})
df = pd.DataFrame({'time': t})

# Generate the data
for a in a_s:
    for s0 in s0_s:
        sol_true = solve_ivp(foo_true, (t[0], t[-1]), (s0,), t_eval=t, args=(a,))
        df_true[f'{a}_{s0}'] = sol_true.y.T
        
for a in a_s:
    for s0 in s0_s:
        sol = solve_ivp(foo, (t[0], t[-1]), (s0,), t_eval=t, args=(a,))    
        df[f'{a}_{s0}'] = sol.y.T

In the above code, df_true is a data frame containing the true dynamics of the system and df the dynamics with a discrepancy.

The figure bellow shows the an example in discrepancy of the data:

Given that I know part of the physics, how can I model $\lambda$?

Does anyone know of a paper/repo with a conceptually similar example that I can a look into?

Best regards

Answer 1

If you know $s(t)$ then you can express the unknown term $\gamma(t) = \lambda(t) + 4a$ as

$$\gamma(t) = 2s(t) + \frac{\text{d}s(t)}{\text{d}t} \approx 2s(t) + \frac{\Delta s(t)}{\Delta t} $$

The approximation works if there is little noise and if the sampling frequency is high (not too much change of the slope within a time period $\Delta t$).

With your example the true model is:

$$s(t) = 2a + (y_0 + 0.125 - 2a) \exp(-2t) + 0.125 \sin(2t) - 0.125 \cos(2t)$$

And an example with $a = 3$ and $y_0 = 3$ looks like

Then if we plot $2s(t) + \frac{\Delta s(t)}{\Delta t}$ we get something that resembles your original sine wave

If you would know the term $a$ then you can subtract it to get $\lambda(t) = \gamma(t) - 4a$. Otherwise you can not know $\gamma(t)$ because you can always subtract or add a constant value to it and you do not know which is the case without any additional information.

---

How you want to estimate/model that final result is a matter that depends on your use case. Note the saying "all models are wrong, but some are useful". You can not sensibly just start modelling with an arbitrary tool without knowing what 'useful' is.

You need context. A question 'how do I model' requires an explanation of the goal why one needs a model in the first place. Why do you need a model if you have the real data without any noise? What is the point of modelling it? With answers to such questions you can follow it up with describing what properties the model should have.

---

Below is some R code to generate the above two figures. You can also make it work with noise, but then the parameters for the filter that is used to compute the derivative need to be changed.

Alternatively, if you know the functional form of $\lambda(t)$ or a family that it belongs to, then you can potentially solve the differential equation and describe $\lambda(t)$ in terms of $s(t)$ without the need for differentiation.

Also, instead of differentiation, one may also use integration (which is less sensitive to noise), but the plotted result $\int_0^t \gamma(x) dx$ is not always easy to interpret. An example of using integration is in the use of fitting a function that is a sum of exponential terms $\sum a_i exp(-b_i t)$ to find a starting condition for other algorithms.

### generate data
dt = 5/1000
t = seq(0, 5, dt)
a = 3
y0 = 3
s = 2*a + (y0 + 0.125 - 2*a) * exp(-2*t) + 0.125 * sin(2*t) - 0.125 * cos(2*t)

### plot data
plot(t,s, type = "l",
     main = expression(s(t) == 2*a + (y[0] + 0.125 - 2*a) * exp(-2*t) + 0.125 * sin(2*t) - 0.125 * cos(2*t)),
     cex.main = 0.8)

### compute lambda
### savgol is a Savitzky Golay filter that we use to compute a derivative
gamma = 2*s + pracma::savgol(s, fl = 3, forder = 1, dorder = 1)/dt

### plotting
range = 2:(length(t)-1) ### trim first and last case where derivative is not computed
plot(t[range], gamma[range], type = "l", ylim = c(11,13),
     main = expression(gamma(t) == 2*s(t) + over(Delta * s(t), Delta * t)),
cex.main = 0.8, ylab = "gamma", xlab = "t")