13

I am trying to implement 2D Laplace inversions through the use of Fourier Series in R. My code worked in the single inversion case but not for double inversion. Any help would be very much appreciated. Thank you!

Below is the code in R for a simple function f(x,t) = exp(-3(x+t)) with Double Laplace Transform 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))}
4

0 回答 0