# Modeling uncertainty from known physics

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, 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, 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

New contributor
nop nop is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct.
• You need to explain what you mean be generating data from this differential equation, what the actual data looks like, and if there is any statistical model involved May 25 at 11:37
• In the post, you have code that will solve the ODE with ground truth dynamic df_true and with some discrepancy df for a range of initial conditions and variations of $a$. I would like to use ML to infer a correction term between the data in df and df_true. May 25 at 11:41
• People shouldn't have to decipher your code in order to understand the question. You should explain more clearly what is the input data and what you are trying to estimate from it. May 25 at 11:53
• It may possibly better if you use a real example instead of a toy example. It is unclear what 'machine learning' is gonna add when you already have the unknown term by simply rewriting the equation where the unknown term is expressed in terms of the known terms. 2 days ago
• The answer to 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 what properties should the model have. [Think about the saying "all models are wrong, but some are useful", you can not sensibly just start modelling with any tool without knowing what 'useful' is. You need context.] 2 days ago

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.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")