# Numerical Laplace inversion with fmincon

3 views (last 30 days)
Neng-Chuan Tien on 7 Nov 2020
Commented: Neng-Chuan Tien on 14 May 2022
Hello
I am using fmincon to fit my experiment data, but now I encountered a problem by using the numerical Laplace m file developed by Tucker McClure.
The m file is:
function ilt = talbot_inversion(f_s, t, M)
% ilt = talbot_inversion(f_s, t, [M])
%
% Returns an approximation to the inverse Laplace transform of function
% handle f_s evaluated at each value in t (1xn) using Talbot's method as
% summarized in the source below.
%
% This implementation is very coarse; use talbot_inversion_sym for better
% precision. Further, please see example_inversions.m for discussion.
%
% f_s: Handle to function of s
% t: Times at which to evaluate the inverse Laplace transformation of f_s
% M: Optional, number of terms to sum for each t (64 is a good guess);
% highly oscillatory functions require higher M, but this can grow
% unstable; see test_talbot.m for an example of stability.
%
te, Joseph, and Ward Whitt. "A Unified Framework for Numerically
% Inverting Laplace Transforms." INFORMS Journal of Computing, vol. 18.4
% (2006): 408-421. Print.
%
% The paper is also online: http://www.columbia.edu/~ww2040/allpapers.html.
%
% Tucker McClure
% Copyright 2012, The MathWorks, Inc.
% Make sure t is n-by-1.
if size(t, 1) == 1
t = t';
elseif size(t, 2) > 1
error('Input times, t, must be a vector.');
end
% Set M to 64 if user didn't specify an M.
if nargin < 3
M = 64;
end
% Vectorized Talbot's algorithm
k = 1:(M-1); % Iteration index
% Calculate delta for every index.
delta = zeros(1, M);
delta(1) = 2*M/5;
delta(2:end) = 2*pi/5 * k .* (cot(pi/M*k)+1i);
% Calculate gamma for every index.
gamma = zeros(1, M);
gamma(1) = 0.5*exp(delta(1));
gamma(2:end) = (1 + 1i*pi/M*k.*(1+cot(pi/M*k).^2)-1i*cot(pi/M*k))...
.* exp(delta(2:end));
% Make a mesh so we can do this entire calculation across all k for all
% given times without a single loop (it's faster this way).
[delta_mesh, t_mesh] = meshgrid(delta, t);
gamma_mesh = meshgrid(gamma, t);
% Finally, calculate the inverse Laplace transform for each given time.
ilt = 0.4./t .* sum(real( gamma_mesh ...
.* arrayfun(f_s, delta_mesh./t_mesh)), 2);
end
and my code is
clear all;
close all;
clc;
Data = ...
[1 0.003322259
2 0.003599114
3 0.003875969
4 0.003322259
5 0.003599114
6 0.003875969
7 0.003322259
8 0.010243632
9 0.127906977
10 0.470376523
11 0.724529347
12 0.63538206
13 0.929678848
14 0.973698782
15 0.978682171
16 0.989202658
17 1.006921373
18 0.996954596
19 0.991971207
20 1.005813953
21 1.029900332
22 1.039590255
23 0.98200443
24 1.003322259
25 0.975636766];
xdata = Data(:,1);
ydata = Data(:,2);
plot(xdata,ydata,'ro')
F = 1.0;
rho = 1.4173;
Kd = 0.0;
V = 5.0;
D = 1.1265;
ZX = 30.0;
theta = 0.475296;
alpha = 0.00001;
beta1 = 1+F*rho*Kd/theta;
beta2 = 1.0;
beta3 = (1-F)*rho*Kd;
x = [V D]
f_s = @(s) exp(x(1)*ZX/2/x(2))/s*exp(-ZX/sqrt(x(2))*sqrt(beta1*s+x(1)^2/4/x(2)+alpha* ...
beta3/theta+alpha^2*beta3/theta/(alpha+beta2*s)))
objective = @(x) sum((talbot_inversion(f_s, xdata) - ydata) .^ 2)
x0 = [3 1];
A = [];
b = [];
Aeq = [];
beq = [];
lb = [];
ub = [];
[x, fval] = fmincon(objective, x0, A, b, Aeq, beq, lb, ub)
'fmincon' does not work since variable x is not include in f_s. Could someone comment to make modifications either in the m-file or in my code that can include variable x in the f_s function or any forward comments? Thank you very much.
Nengchuan
Lateef Adewale Kareem on 29 Apr 2022
fun = @(s, x)exp(x(1)*ZX/2/x(2))/s*exp(-ZX/sqrt(x(2))*sqrt(beta1*s+x(1)^2/4/x(2)+alpha* ...
beta3/theta+alpha^2*beta3/theta/(alpha+beta2*s)))
objective = @(x) sum((talbot_inversion(@(s)fun(s,x), xdata) - ydata) .^ 2)

Lateef Adewale Kareem on 29 Apr 2022
I am sorry, I didnt see the problem early. This is the complete code
clear all;
close all;
clc;
Data = ...
[1 0.003322259
2 0.003599114
3 0.003875969
4 0.003322259
5 0.003599114
6 0.003875969
7 0.003322259
8 0.010243632
9 0.127906977
10 0.470376523
11 0.724529347
12 0.63538206
13 0.929678848
14 0.973698782
15 0.978682171
16 0.989202658
17 1.006921373
18 0.996954596
19 0.991971207
20 1.005813953
21 1.029900332
22 1.039590255
23 0.98200443
24 1.003322259
25 0.975636766];
xdata = Data(:,1); ydata = Data(:,2);
F = 1.0;
rho = 1.4173;
Kd = 0.0;
V = 5.0;
D = 1.1265;
ZX = 30.0;
theta = 0.475296;
alpha = 0.00001;
beta1 = 1+F*rho*Kd/theta;
beta2 = 1.0;
beta3 = (1-F)*rho*Kd;
fun = @(s, x)exp(x(1)*ZX/2/x(2))/s*exp(-ZX/sqrt(x(2))*sqrt(beta1*s+ ...
x(1)^2/4/x(2)+alpha*beta3/theta+alpha^2*beta3/theta/(alpha+beta2*s)));
objective = @(x) sum((talbot_inversion(@(s)fun(s,x), xdata) - ydata) .^ 2);
x0 = [V, D];
A = [];
b = [];
Aeq = [];
beq = [];
lb = [];
ub = [];
[x, fval] = fmincon(objective, x0, A, b, Aeq, beq, lb, ub);
Local minimum possible. Constraints satisfied. fmincon stopped because the size of the current step is less than the value of the step size tolerance and constraints are satisfied to within the value of the constraint tolerance.
ymodel_initialguess = talbot_inversion(@(s)fun(s,x0), xdata);
ymodel = talbot_inversion(@(s)fun(s,x), xdata);
plot(xdata,ydata,'ko'); hold on
plot(xdata, ymodel_initialguess, 'linewidth', 2);
plot(xdata, ymodel, 'linewidth', 2);
xlabel('xdata', 'interpreter', 'Latex');
ylabel('ydata', 'interpreter', 'Latex');
set(gca, 'FontSize', 15, 'TickLabelInterpreter', 'latex');
legend('Experiment', 'Model InitialGuess', 'FinalModel', ...
'interpreter', 'Latex', 'location', 'best'); function ilt = talbot_inversion(f_s, t, M)
% ilt = talbot_inversion(f_s, t, [M])
%
% Returns an approximation to the inverse Laplace transform of function
% handle f_s evaluated at each value in t (1xn) using Talbot's method as
% summarized in the source below.
%
% This implementation is very coarse; use talbot_inversion_sym for better
% precision. Further, please see example_inversions.m for discussion.
%
% f_s: Handle to function of s
% t: Times at which to evaluate the inverse Laplace transformation of f_s
% M: Optional, number of terms to sum for each t (64 is a good guess);
% highly oscillatory functions require higher M, but this can grow
% unstable; see test_talbot.m for an example of stability.
% Joseph, and Ward Whitt. "A Unified Framework for Numerically
% Inverting Laplace Transforms." INFORMS Journal of Computing, vol. 18.4
% (2006): 408-421. Print.
%
% The paper is also online: http://www.columbia.edu/~ww2040/allpapers.html.
%
% Tucker McClure
% Copyright 2012, The MathWorks, Inc.
% Make sure t is n-by-1.
if size(t, 1) == 1
t = t';
elseif size(t, 2) > 1
error('Input times, t, must be a vector.');
end
% Set M to 64 if user didn't specify an M.
if nargin < 3
M = 64;
end
% Vectorized Talbot's algorithm
k = 1:(M-1); % Iteration index
% Calculate delta for every index.
delta = zeros(1, M);
delta(1) = 2*M/5;
delta(2:end) = 2*pi/5 * k .* (cot(pi/M*k)+1i);
% Calculate gamma for every index.
gamma = zeros(1, M);
gamma(1) = 0.5*exp(delta(1));
gamma(2:end) = (1 + 1i*pi/M*k.*(1+cot(pi/M*k).^2)-1i*cot(pi/M*k))...
.* exp(delta(2:end));
% Make a mesh so we can do this entire calculation across all k for all
% given times without a single loop (it's faster this way).
[delta_mesh, t_mesh] = meshgrid(delta, t);
gamma_mesh = meshgrid(gamma, t);
% Finally, calculate the inverse Laplace transform for each given time.
ilt = 0.4./t .* sum(real( gamma_mesh ...
.* arrayfun(f_s, delta_mesh./t_mesh)), 2);
end
Neng-Chuan Tien on 14 May 2022
Thank you very much!!