Двойная инверсия Лапласа в R с серией Фурье - программирование

Двойная инверсия Лапласа в R с серией Фурье

Я пытаюсь реализовать 2D-образы Лапласа с помощью рядов Фурье в R. Мой код работал в случае единственной инверсии, но не для двойной инверсии. Любая помощь будет очень высоко ценится. Спасибо!

Ниже приведен код в R для простой функции f(x,t) = exp(-3(x+t)) с двойным преобразованием Лапласа F(p,s) = 1/((p+3)*(s+3)).

# General implementation of 2D Laplace Inversion

# How to use this?
# 1. Change the Laplace Transform function (F)
# 2. Change the Original function(f_1)

# Double Laplace Transform function, parameters p and s
F <- function(p,s) {1/((p+3)*(s+3))}

# Adjust parameters here
a <- 0.1
n <- 1000
# Note that the approximation works accurately up to x=T/2 only
T <- 10

# First inversion - p
    F_1 <- function(x,s,k) {Re(F(complex(real=a, imaginary=k*pi/T),s))*cos(k*pi*x/T)}
    # Summing over k from 1 to n
    F_2 <- function(x,s) {sapply(1:n, F_1, x=x, s=s)}
    # Approximate Function using Fourier Series 
    G <- function(x,s) {2*exp(a*x)/T*(0.5*Re(F(a,s))+sum(F_2(x,s)))}

# Second inversion - s
    G_1 <- function(x,t,k) {Re(G(x,complex(real=a, imaginary=k*pi/T)))*cos(k*pi*t/T)}
    # Summing over l from 1 to n
    G_2 <- function(x,t) {sapply(1:n, G_1, x=x, t=t)}
    # Approximate Function using Fourier Series 
    f2 <- function(x,t) {2*exp(a*t)/T*(0.5*Re(G(x,a))+sum(G_2(x,t)))}

# The other way round.
# First inversion - s
    F_1 <- function(t,p,k) {Re(F(p, complex(real=a, imaginary=k*pi/T)))*cos(k*pi*t/T)}
    # Summing over k from 1 to n
    F_2 <- function(t,p) {sapply(1:n, F_1, t=t, p=p)}
    # Approximate Function using Fourier Series 
    G <- function(t,p) {2*exp(a*t)/T*(0.5*Re(F(p,a))+sum(F_2(t,p)))}

# Second inversion - p
    G_1 <- function(x,t,k) {Re(G(t,complex(real=a, imaginary=k*pi/T)))*cos(k*pi*x/T)}
    # Summing over l from 1 to n
    G_2 <- function(x,t) {sapply(1:n, G_1, x=x, t=t)}
    # Approximate Function using Fourier Series 
    f2 <- function(x,t) {2*exp(a*x)/T*(0.5*Re(G(a,t))+sum(G_2(x,t)))}

# Original Function
f1 <- function(x,t) {exp(-3*(x+t))}
4b9b3361