Skip to content

come valutare f su valori X "fittizi"? #3

@AndreaPi

Description

@AndreaPi

Quasi tutti i metodi numerici più accurati di Eulero esplicito richiedono di valutare dynamic_f (il RHS del sistema SIR) per dei valori X "fittizi" che non corrispondono ad alcun valore reale della variabile X. Per esempio, metodo di Heun, eq. 11.10:

http://www.ceri.memphis.edu/people/echoi2/ceri8315/Quarteroni-ComputMath-Ch11.pdf

Usando come template il codice di euler,
https://github.com/enricomeloni/covid-tools/blob/baeb5ad140a9742ebab466fbf497d8c6aaccdd70/torch_euler.py#L36

otterremmo per il metodo di Heun:

def Heun (f, omega, time_grid):
    """

    :param f: function describing the differential equations
    :param omega: function returning values from -inf to 0 as a N-Dim tuple
    :param time_grid: 1-Dim tensor representing the time-grid on which to integrate
    :return: 1-Dim tensor the same size as time_grid with values computed on the time grid

    NOTE: only expected to be second-order accurate if dt is constant
    """

    y0 = torch.tensor([omega(0)])
    time_grid = time_grid.to(y0[0])
    values = y0


    for i in range(0, time_grid.shape[0] - 1):
        t_i = time_grid[i]
        t_next = time_grid[i+1]
        y_i = values[i]
        dt = torch.tensor([t_next - t_i])
        t_series = TimeSeries(time_grid, values, omega)
        f1 = torch.stack(tuple(f_t.unsqueeze(0) for f_t in f(t_i, t_series, dt)), dim=1)
        f2 = torch.stack(tuple(f_t.unsqueeze(0) for f_t in f(t_next, ...., dt)), dim=1)
        dy = 0.5*dt*(f1+f2)
        y_next = y_i + dy
        #y_next = y_next.unsqueeze(0)
        values = torch.cat((values, y_next), dim=0)

    return values

Nel calcolo di f2, ho bisogno di bisogno di valutare f nel punto (t_next, y_i + f1)...ma dato che y_i + f1 è un valore che non viene mai realmente assunto dalla funzione, non c'è in t_series. Come faccio a passarlo?

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions