Reverse problem of finding time-varying parameters of an ODE with the help of solution data.

Question

Anshuman on 23 Nov 2022

0
Link

Direct link to this question

https://www.mathworks.com/matlabcentral/answers/1860138-reverse-problem-of-finding-time-varying-parameters-of-an-ode-with-the-help-of-solution-data

Commented: Anshuman on 23 Nov 2022

Any way in which I can find the time-varying parameters of an ODE ; The solution to which is know in the form of time-series data.

Background:

Example: SIR model ode: as following.

% The SIR model differential equations.

def deriv(y, t, beta, gamma):

S, I, R = y

dSdt = -beta * S * I

dIdt = beta * S * I - gamma * I

dRdt = gamma * I

return dSdt, dIdt, dRdt

Then, I use odeint(deriv, y0, t, args=( beta, gamma)) (In python) to solve for S, I and R using minimisation of these parameters.

Question

But I want to ask, is there any way in which a reverse problem can be constructed such that; I have the data for S, I and R as time-sereis. and then we can calculate 'beta' and 'gamma' paramter from there?

Note beta and gamma should be time dependent.They can be assumed as sum of natural cubic splines; Where N_k(t), k=1 to K, are K natural cubic spline basis functions evaluated at K-2 equally spaced knots in addition to the boundary knots .

Can this problem be solved on matlab?

Thanks.

7 Comments
Show 5 older commentsHide 5 older comments

Anshuman on 23 Nov 2022

Open in MATLAB Online

I have written a code for estimating 'beta' and ;delta' values for the whole time step; that is it will return a constant beta and delta value for the whole timespan in python.

Code:

(a) Define the ODE model of SIR with cosntant parameter 'beta' and 'gamma'

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

I0, R0 = 0.10798003802281400, 0.0020386203150461700
S0 = 1 - I0 - R0
t = np.linspace(0, 146, 146)

# The SIR model differential equations.

def deriv(y, t, beta, gamma):
    S, I, R = y
    dSdt = -beta * S * I 
    dIdt = beta * S * I  - gamma * I
    dRdt = gamma * I
    return dSdt, dIdt, dRdt

I_d = x[:,1] # from dataset, we will finally try to fit the SIR plot with this data 

(b) Define fn(x) for minimising the loss with respect to the data of 'I'

def fn(x):
# parameters unwrapped
    beta, gamma = x;   
    I0, R0 = 0.10798003802281400, 0.0020386203150461700
    S0 = 1 - I0 - R0
    y0 = S0, I0, R0
    ret = odeint(deriv, y0, t_d, args=(beta, gamma))
    S, I, R = ret.T
    error = sum((I-I_d)*(I-I_d))
    return error

# initialise with current best guess
init_x = [0.1, 0.05]
res = minimize(fn, init_x, method='Nelder-Mead', tol=1e-8)
fn(res.x) # calculate final error
res.x # calculate final parameters