Numerical Laplace inversion with fmincon

2 次查看（过去 30 天）
Neng-Chuan Tien
Neng-Chuan Tien 2020-11-7
评论： Neng-Chuan Tien 2022-5-14
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
  1 个评论
Lateef Adewale Kareem
Lateef Adewale Kareem 2022-4-29
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
Ran in:
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

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Systems of Nonlinear Equations 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by