Problem with system of ODEs

[sebastian-griego]

I'm having trouble with getting the right solution in a forward problem with a system of ODEs. I have tried to follow the examples, but I do not get the right solution regardless. This is similar to #1719. I would appreciate any help.

import re
import numpy as np
import matplotlib.pyplot as plt
import deepxde as dde
from deepxde.backend import tf
from scipy.integrate import solve_ivp, odeint
import math
import pandas as pd
from io import StringIO
from scipy.interpolate import interp1d

# Initial Conditions
T0 = 100000
I0 = 0
L0 = 0
V0 = 200

# Parameters
d = 0.01
beta = 0.0000002
f = 0.001
a = 0.2
delta_I = 1 
delta_L = 0.0039
c = 23
p = 4000
lambda_param = 1000  # T0 * d
Omega = 1

t_initial = 0
t_final = 100

geom = dde.geometry.TimeDomain(t_initial, t_final)

Y0 = [T0, I0, L0, V0]

def system(Y, t):
    T, I, L, V = Y
    dT_dt = lambda_param - d * T - Omega * beta * V * T
    dI_dt = (1 - f) * Omega * beta * V * T + a * L - delta_I *   I
    dL_dt = f * Omega * beta * V * T - a * L - delta_L * L
    dV_dt = Omega * p * I - c * V
    return [dT_dt, dI_dt, dL_dt, dV_dt]

t = np.linspace(t_initial, t_final, 1000)

solution = odeint(system, Y0, t)
T, I, L, V = solution.T

interp_T = interp1d(t, T, kind='cubic', fill_value='extrapolate')
interp_I = interp1d(t, I, kind='cubic', fill_value='extrapolate')
interp_L = interp1d(t, L, kind='cubic', fill_value='extrapolate')
interp_V = interp1d(t, V, kind='cubic', fill_value='extrapolate')

def func(t):
    T_val = interp_T(t)
    I_val = interp_I(t)
    L_val = interp_L(t)
    V_val = interp_V(t)
    return np.hstack((T_val, I_val, L_val, V_val))

def ode(t, Y):
    T = Y[:, 0:1]
    I = Y[:, 1:2]
    L = Y[:, 2:3]
    V = Y[:, 3:4]
    
    dT_dt = dde.grad.jacobian(Y, t, i=0)
    dI_dt = dde.grad.jacobian(Y, t, i=1)
    dL_dt = dde.grad.jacobian(Y, t, i=2)
    dV_dt = dde.grad.jacobian(Y, t, i=3)
    
    return [
        dT_dt - (lambda_param - d * T - Omega * beta * V * T),
        dI_dt - ((1 - f) * Omega * beta * V * T + a * L - delta_I * I),
        dL_dt - (f * Omega * beta * V * T - a * L - delta_L * L),
        dV_dt - (Omega * p * I - c * V)
    ]

def boundary(_, on_initial):
    return on_initial

ic_T = dde.icbc.IC(geom, lambda x: T0, boundary, component=0)
ic_I = dde.icbc.IC(geom, lambda x: I0, boundary, component=1)
ic_L = dde.icbc.IC(geom, lambda x: L0, boundary, component=2)
ic_V = dde.icbc.IC(geom, lambda x: V0, boundary, component=3)

data = dde.data.PDE(
    geom,
    ode,
    [ic_T, ic_I, ic_L, ic_V],
    num_domain=3000,
    num_boundary=2,
    solution=func,
    num_test=1000
)

neurons = 64
layers = 3
activation = "relu"
initializer = "Glorot normal"
net = dde.nn.FNN([1] + [neurons] * layers + [4], activation, initializer)

model = dde.Model(data, net)

model.compile("adam", lr=0.001, loss_weights=[1e-6, 1, 1, 1e-8, 1e-10, 1, 1, 1e-4])
losshistory, train_state = model.train(iterations=25000, display_every=5000)