f =
plot after taking inversion Laplace to val function and then divide it on val0 function for steady case, it should give 1 as t tends to infty
Show older comments
plot after taking inversion Laplace to val function and then divide it on val0 function for steady case, it should give 1 as t tends to infty
function main_calculation()
clc; close all;
t_values = logspace(0, 2, 100);
beta1_vals = [0.001, 1, 3, 50]; % 4 specific curves to compare
params = struct('c',0.0000001, 's1',0.00001, 's2',0.1, 'rho1',0, 'beta2',10000, 'b',2,'j',0.1);
figure; hold on;
colors = ['r', 'g', 'b', 'm']; % Three colors for the three curves
for k = 1:length(beta1_vals)
params.beta1 = beta1_vals(k);
u_s1 = zeros(size(t_values));
u_steady = zeros(size(t_values));
for i = 1:length(t_values)
% Compute both inversions at the current beta1
u_s1(i) = talbot_inversion(@(s) U_func_s1(s, params), t_values(i));
u_steady(i) = talbot_inversion(@(s) U_func_steady(s, params), t_values(i));
end
% Calculate ratio
ratio = u_s1 ./ (u_steady + eps);
% Plot
plot(t_values, ratio, 'Color', colors(k), 'LineWidth', 2, ...
'DisplayName', sprintf('\\beta_1=%.2f', beta1_vals(k)));
end
set(gca, 'XScale', 'log');
xlabel('Time (t)'); ylabel('Ratio (U_{osc}/U_{std})');
legend('Location', 'best'); grid on;
% title('Effect of Micropolar Coupling on Dispersion Ratio');
end
% --- Rename your existing logic to avoid conflicts ---
function val = U_func_s1(sigma, p)
% Copy your logic from the original calc_s1.m U_function here
% Ensure it returns 'val'
% Unpack
b = p.b; c = p.c; s1 = p.s1; s2 = p.s2;
beta1 = p.beta1; beta2 = p.beta2; j = p.j; rho1 = p.rho1;
% Define S (sum) and P (product)
S = (4*c*s1 + (1 + j*s1*(c + 1))*sigma) / (c + 1);
P = (sigma*s1*(4*c + j*sigma)) / (c + 1);
% Calculate alpha1 and alpha2
alpha1 = sqrt((S - sqrt(S^2 - 4*P)) / 2);
alpha2 = sqrt((S + sqrt(S^2 - 4*P)) / 2);
%A1,B1,C1,D1,E1,F1
a11= 2; a12= 2 ;
a13= 2*besselk(0.3e1 / 0.2e1, alpha1);a14= 2 * besselk(0.3e1 / 0.2e1, alpha2) ;
a15= 2*besseli(0.3e1 / 0.2e1, alpha1);a16= 2 * besseli(0.3e1 / 0.2e1, alpha2);
a21= -(beta1 + 4 * c + 2) / beta1 ;
a22= 2 * (beta1 - 2 * c + 2) / beta1 ;
a23=(-alpha1 * (beta1 + 2) / beta1 * besselk(0.1e1 / 0.2e1, alpha1) - (beta1 + 4 * c + sigma + 2) / beta1 * besselk(0.3e1 / 0.2e1, alpha1));
a24=(-alpha2 * (beta1 + 2) / beta1 * besselk(0.1e1 / 0.2e1, alpha2) - (beta1 + 4 * c + sigma + 2) / beta1 * besselk(0.3e1 / 0.2e1, alpha2));
a25=(alpha1 * (beta1 + 2) / beta1 * besseli(0.1e1 / 0.2e1, alpha1) - (beta1 + 4 * c + sigma + 2) / beta1 * besseli(0.3e1 / 0.2e1, alpha1)) ;
a26=(alpha2 * (beta1 + 2) / beta1 * besseli(0.1e1 / 0.2e1, alpha2) - (beta1 + 4 * c + sigma + 2) / beta1 * besseli(0.3e1 / 0.2e1, alpha2));
%a31=0;a32=0;
a33= ((alpha1 ^ 2 * beta2 * c * s1 * s2 + alpha1 ^ 2 * beta2 * s1 * s2 - alpha1 ^ 2 * c * s1 - beta2 * s1 * s2 * sigma + alpha1 ^ 2 * c - alpha1 ^ 2 * s1 + alpha1 ^ 2 + s1 * sigma - sigma) / beta2 / s1 / c / s2 * besselk(0.3e1 / 0.2e1, alpha1) / 0.2e1 + 0.1e1 / beta2 / s2 / c * (alpha1 ^ 2 * c + alpha1 ^ 2 - sigma) * besselk(0.5e1 / 0.2e1, alpha1) * alpha1 / 0.2e1) ;
a34= ((alpha2 ^ 2 * beta2 * c * s1 * s2 + alpha2 ^ 2 * beta2 * s1 * s2 - alpha2 ^ 2 * c * s1 - beta2 * s1 * s2 * sigma + alpha2 ^ 2 * c - alpha2 ^ 2 * s1 + alpha2 ^ 2 + s1 * sigma - sigma) / beta2 / s1 / c / s2 * besselk(0.3e1 / 0.2e1, alpha2) / 0.2e1 + 0.1e1 / beta2 / s2 / c * (alpha2 ^ 2 * c + alpha2 ^ 2 - sigma) * besselk(0.5e1 / 0.2e1, alpha2) * alpha2 / 0.2e1) ;
a35= ((alpha1 ^ 2 * beta2 * c * s1 * s2 + alpha1 ^ 2 * beta2 * s1 * s2 - alpha1 ^ 2 * c * s1 - beta2 * s1 * s2 * sigma + alpha1 ^ 2 * c - alpha1 ^ 2 * s1 + alpha1 ^ 2 + s1 * sigma - sigma) / beta2 / s1 / c / s2 * besseli(0.3e1 / 0.2e1, alpha1) / 0.2e1 - 0.1e1 / beta2 / s2 / c * (alpha1 ^ 2 * c + alpha1 ^ 2 - sigma) * alpha1 * besseli(0.5e1 / 0.2e1, alpha1) / 0.2e1) ;
a36= ((alpha2 ^ 2 * beta2 * c * s1 * s2 + alpha2 ^ 2 * beta2 * s1 * s2 - alpha2 ^ 2 * c * s1 - beta2 * s1 * s2 * sigma + alpha2 ^ 2 * c - alpha2 ^ 2 * s1 + alpha2 ^ 2 + s1 * sigma - sigma) / beta2 / s1 / c / s2 * besseli(0.3e1 / 0.2e1, alpha2) / 0.2e1 - 0.1e1 / beta2 / s2 / c * (alpha2 ^ 2 * c + alpha2 ^ 2 - sigma) * besseli(0.5e1 / 0.2e1, alpha2) * alpha2 / 0.2e1) ;
a41=2 / b ^ 3; a42= 2;
a43= 2 * b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha1) ;
a44= 2 * b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha2);
a45= 2 * b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha1) ;
a46= 2 * b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha2);
a51= -1 / b ^ 3 ;a52= 2 ;
a53= (-b ^ (-0.1e1 / 0.2e1) * alpha1 * besselk(0.1e1 / 0.2e1, b * alpha1) - b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha1)) ;
a54= (-b ^ (-0.1e1 / 0.2e1) * alpha2 * besselk(0.1e1 / 0.2e1, b * alpha2) - b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha2)) ;
a55= (b ^ (-0.1e1 / 0.2e1) * alpha1 * besseli(0.1e1 / 0.2e1, b * alpha1) - b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha1)) ;
a56= (b ^ (-0.1e1 / 0.2e1) * alpha2 * besseli(0.1e1 / 0.2e1, b * alpha2) - b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha2)) ;
%a61=0; a62=0;
a63= 0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha1 ^ 2 * b ^ 2 * c + alpha1 ^ 2 * b ^ 2 - b ^ 2 * sigma) * besselk(0.3e1 / 0.2e1, b * alpha1)/ 0.2e1;
a64= 0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha2 ^ 2 * b ^ 2 * c + alpha2 ^ 2 * b ^ 2 - b ^ 2 * sigma) * besselk(0.3e1 / 0.2e1, b * alpha2) / 0.2e1;
a65= 0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha1 ^ 2 * b ^ 2 * c + alpha1 ^ 2 * b ^ 2 - b ^ 2 * sigma) * besseli(0.3e1 / 0.2e1, b * alpha1) / 0.2e1 ;
a66= 0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha2 ^ 2 * b ^ 2 * c + alpha2 ^ 2 * b ^ 2 - b ^ 2 * sigma) * besseli(0.3e1 / 0.2e1, b * alpha2) / 0.2e1;
% Construct the 6x6 matrix A and column vector B
A = [a11, a12, a13, a14, a15, a16;
a21, a22, a23, a24, a25, a26;
0, 0, a33, a34, a35, a36;
a41, a42, a43, a44, a45, a46;
a51, a52, a53, a54, a55, a56;
0, 0, a63, a64, a65, a66];
B = [1; -1; 0; 0; 0; 0]; % Example: solve Ax=B for unit source
x = A \ B;
x1 = x(1);
val =( 3/(4*pi*sigma^2*(rho1 - 1 - 3*x1)));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function val0 = U_func_steady(sigma, p)
% Copy your logic from the original calc_steady.m U_function here
% Ensure it returns 'val0'
% Unpack
b = p.b; c = p.c; s1 = p.s1; s2 = p.s2;
beta1 = p.beta1; beta2 = p.beta2;
% Define S (sum) and P (product)
alpha = sqrt(4*c*s1 /(c + 1));
a11= 2; a12= 2 ; a13= 2 ; a14=2; a15=2 * besselk(0.3e1 / 0.2e1, alpha); a16=2 * besseli(0.3e1 / 0.2e1, alpha);
a21 = -(beta1 + 4 * c + 2) / beta1 ;
a22=2 * (beta1 - 2 * c + 2) / beta1 ;
a23=2 * (2 * beta1 - 7 * c - 1) / beta1;
a24=(beta1 - 2 * c + 4) / beta1 ;
a25=(-alpha * (beta1 + 2) / beta1 * besselk(0.1e1 / 0.2e1, alpha) - (beta1 + 4 * c + sigma + 2) / beta1 * besselk(0.3e1 / 0.2e1, alpha));
a26=(alpha * (beta1 + 2) / beta1 * besseli(0.1e1 / 0.2e1, alpha) - (beta1 + 4 * c + sigma + 2) / beta1 * besseli(0.3e1 / 0.2e1, alpha));
%a31=0;a32=0;a33=0;a34=0;
a35= ((alpha ^ 2 * beta2 * c * s1 * s2 + alpha ^ 2 * beta2 * s1 * s2 - alpha ^ 2 * c * s1 - beta2 * s1 * s2 * sigma + alpha ^ 2 * c - alpha ^ 2 * s1 + alpha ^ 2 + s1 * sigma - sigma) / beta2 / s1 / c / s2 * besselk(0.3e1 / 0.2e1, alpha) / 0.2e1 + 0.1e1 / beta2 / s2 / c * (alpha ^ 2 * c + alpha ^ 2 - sigma) * besselk(0.5e1 / 0.2e1, alpha) * alpha / 0.2e1) ;
a36= ((alpha ^ 2 * beta2 * c * s1 * s2 + alpha ^ 2 * beta2 * s1 * s2 - alpha ^ 2 * c * s1 - beta2 * s1 * s2 * sigma + alpha ^ 2 * c - alpha ^ 2 * s1 + alpha ^ 2 + s1 * sigma - sigma) / beta2 / s1 / c / s2 * besseli(0.3e1 / 0.2e1, alpha) / 0.2e1 - 0.1e1 / beta2 / s2 / c * (alpha ^ 2 * c + alpha ^ 2 - sigma) * besseli(0.5e1 / 0.2e1, alpha) * alpha / 0.2e1) ;
a41=2 / b ^ 3; a42=2 ; a43=2 * b ^ 2 ; a44=2 / b ;
a45=2 * b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha);
a46=2 * b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha);
a51= -1 / b ^ 3; a52= 2 ; a53=4 * b ^ 2 ;a54=1 / b ;
a55=(-b ^ (-0.1e1 / 0.2e1) * alpha * besselk(0.1e1 / 0.2e1, b * alpha) - b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha)) ;
a56=(b ^ (-0.1e1 / 0.2e1) * alpha * besseli(0.1e1 / 0.2e1, b * alpha) - b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha)) ;
%a61=0;a62=0;a63=0;a64=0;
a65= 0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha ^ 2 * b ^ 2 * c + alpha ^ 2 * b ^ 2 - b ^ 2 * sigma) * besselk(0.3e1 / 0.2e1, b * alpha) / 0.2e1;
a66= 0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha ^ 2 * b ^ 2 * c + alpha ^ 2 * b ^ 2 - b ^ 2 * sigma) * besseli(0.3e1 / 0.2e1, b * alpha) / 0.2e1;
% Construct the 6x6 matrix A and column vector B
A0 = [a11, a12, a13, a14, a15, a16;
a21, a22, a23, a24, a25, a26;
0, 0, 0, 0, a35, a36;
a41, a42, a43, a44, a45, a46;
a51, a52, a53, a54, a55, a56;
0, 0, 0, 0, a65, a66];
B0 = [1; -1; 0; 0; 0; 0]; % Example: solve Ax=B for unit source
x0 = A0 \ B0;
x0 = x0(1);
% val0 =(-3/(4*pi*(1 + 3*x0)));
val0 =-1/(4*pi*(1 +c)*x0);
end
1 Comment
Why does u_steady depend on t ? E.g. if your inverse Laplace transform is 5 - 5*exp(-2*t), u_s1(t) = 5 - 5*exp(-2*t) and u_steady = 5 (which does not depend on t). Or do I misunderstand your code ?
syms s
f = 10/(s*(s+2))
ilaplace(f)
Answers (1)
ADDENDUM after sidebar conversation --
...
% calculate SS result as functional value at t==inf
u_steady=talbot_inversion(@(s) U_func_s1(s, params), inf);
% Calculate ratio
ratio = u_s1/u_steady;
...
If the function doesn't like inf then
...
% calculate SS result as functional value at t=realmax
u_steady=talbot_inversion(@(s) U_func_s1(s, params), realmax);
% Calculate ratio
ratio = u_s1/u_steady;
...
will be at the largest value that MATLAB double will hold which is the best can approximate infinity if the functional can't handle the IEEE inf.
main_calculation
function main_calculation()
...CODE ELIDED FOR POST BREVITY -- DPB
end
Missing the talbot_inversion function so can't be run as is.
Without digging more deeply to see what the system would reduce to vs time, the empirical way would be to simply use a very large time value as a single point input into the same function for both cases since end steady-state value is a constant.
...
%u_steady = zeros(size(t_values));
for i = 1:length(t_values)
% Compute inversion at the current beta1 and t
u_s1(i) = talbot_inversion(@(s) U_func_s1(s, params), t_values(i));
end
% Compute inversion at the current beta1 and very large t --> inf
u_steady= talbot_inversion(@(s) U_func_s1(s, params), 100*t_values(end));
% Calculate ratio
ratio = u_s1/u_steady; % don't need ./, u_steady is a constant
...
Then you wouldn't need to rename the original and have a second evaluation function.
The above uses an aribrary 100 times the time vector ending value; <as showed previously>, creating variables for the purpose as did for tEnd before to change the time over which the plots were generated would again be the better way; the above has reverted to the fixed 0,10 in the linspace() call.
31 Comments
Shreen El-Sapa
on 8 Jul 2026 at 12:53
Edited: Walter Roberson
on 8 Jul 2026 at 15:27
dpb
on 8 Jul 2026 at 14:05
Torsten
on 8 Jul 2026 at 14:10
For your suggestion, "ratio" will always converge to 1 if u_s1(t) has an asymptotic value as t -> Inf. I don't think that's what @Shreen El-Sapa intended.
I guess I don't follow -- your example above of a specific transfer function
u_s1=@(t)5 - 5*exp(-2*t)
with an asymptotic value of 5 can be generated with
t_values=linspace(0,10,100);
u_steady=u_s1(100*t_values(end))
It isn't at all clear to me without a lot more digging than have time (and patience) to do, to figure out an analytic expression for the subject transfer function here, so just doing it empirically.
What's different from what @Shreen El-Sapa asked for except it isn't analytic? And, certainly the asymptotic value shouldn't be dependent upon t, unless trying to calculate/estimate it from the time solution.
ADDENDUM Maybe you overlooked the "100*" factor on the ending t_value so thought it was the same as the end value calculated, Torsten? That wouldn't be correct and would always end up at unity, yes, but that isn't what was suggested.
Hmmm...
u_steady=u_s1(inf)
will work with the analytic but may not succeed with the OP's function -- I suppose instead of the above arbitrary "100*" trick, one could try realmax instead.
u_steady=u_s1(realmax)
I'll make an addendum to the Answer to note the "enhancement" to get away from the totally empirical multiplication factor and @Shreen El-Sapa can try it out to see if the functional blows up or not.
Shreen El-Sapa
on 8 Jul 2026 at 14:33
Edited: Shreen El-Sapa
on 8 Jul 2026 at 14:35
dpb
on 8 Jul 2026 at 14:52
Yes, that's what my suggestion will do with a calculated value for the asymptotic value from the base function evaluated at essentially t==inf. As noted in the sidebar conversation with @Torsten, using either inf or realmax in place of the empirical multiplier on the ending time would remove that empiricism and give the best estimate you can get from the functional.
Of course, if you do have an analytic expression for the function, then evaluating it is the obvious "correct" way; I presume the whole point of the exercise is that it is only known numerically in which case the suggested solution should work well.
I tried to test the prior suggestion -- doesn't work well with the given functional -- for t=inf, it returns 0 so the ratio --> inf; for t=realmax the functional returned inf so the ratio --> 0
I decided then to explore why the functional didn't behave with large t (first time noticed it was log time, not linear) so expanded the range for the last case (the only one with any major change in values over time) and first the value has a minimum around 50 and then it behaves very badly for large t.
Note the plot is the u value for the last beta_1, not a ratio -- just left the ylabel unchanged.
But, whatever is going on in the functional is not good; it generates condition warnings on any input time, not just large values. Would need to figure out how to fix that before anything else I suspect and then at the present there really isn't an asymptotic value to pick.
warning('off')
main_calculation
function main_calculation()
clc; close all;
t_values = logspace(0, 6, 100);
beta1_vals = [0.001, 1, 3, 50]; % 4 specific curves to compare
params = struct('c',0.0000001, 's1',0.00001, 's2',0.1, 'rho1',0, 'beta2',10000, 'b',2,'j',0.1);
figure; hold on;
colors = ['r', 'g', 'b', 'm']; % Three colors for the three curves
for k = numel(beta1_vals) %1:length(beta1_vals)
params.beta1 = beta1_vals(k);
u_s1 = zeros(size(t_values));
%u_steady = zeros(size(t_values));
for i = 1:length(t_values)
% Compute both inversions at the current beta1
u_s1(i) = talbot_inversion(@(s) U_func_s1(s, params), t_values(i));
end
u_steady= talbot_inversion(@(s) U_func_s1(s, params), 10*t_values(end));
%u_s1
u_steady
% Calculate ratio
ratio = u_s1/u_steady;
% Plot
plot(t_values, u_s1, 'Color', colors(k), 'LineWidth', 2, ...
'DisplayName', sprintf('\\beta_1=%.2f', beta1_vals(k)));
%plot(t_values, ratio, 'Color', colors(k), 'LineWidth', 2, ...
% 'DisplayName', sprintf('\\beta_1=%.2f', beta1_vals(k)));
end
set(gca, 'XScale', 'log');
xlabel('Time (t)'); ylabel('Ratio (U_{osc}/U_{std})');
legend('Location', 'best'); grid on;
% title('Effect of Micropolar Coupling on Dispersion Ratio');
end
% --- Rename your existing logic to avoid conflicts ---
function val = U_func_s1(sigma, p)
% Copy your logic from the original calc_s1.m U_function here
% Ensure it returns 'val'
% Unpack
b = p.b; c = p.c; s1 = p.s1; s2 = p.s2;
beta1 = p.beta1; beta2 = p.beta2; j = p.j; rho1 = p.rho1;
% Define S (sum) and P (product)
S = (4*c*s1 + (1 + j*s1*(c + 1))*sigma) / (c + 1);
P = (sigma*s1*(4*c + j*sigma)) / (c + 1);
% Calculate alpha1 and alpha2
alpha1 = sqrt((S - sqrt(S^2 - 4*P)) / 2);
alpha2 = sqrt((S + sqrt(S^2 - 4*P)) / 2);
%A1,B1,C1,D1,E1,F1
a11= 2; a12= 2 ;
a13= 2*besselk(0.3e1 / 0.2e1, alpha1);a14= 2 * besselk(0.3e1 / 0.2e1, alpha2) ;
a15= 2*besseli(0.3e1 / 0.2e1, alpha1);a16= 2 * besseli(0.3e1 / 0.2e1, alpha2);
a21= -(beta1 + 4 * c + 2) / beta1 ;
a22= 2 * (beta1 - 2 * c + 2) / beta1 ;
a23=(-alpha1 * (beta1 + 2) / beta1 * besselk(0.1e1 / 0.2e1, alpha1) - (beta1 + 4 * c + sigma + 2) / beta1 * besselk(0.3e1 / 0.2e1, alpha1));
a24=(-alpha2 * (beta1 + 2) / beta1 * besselk(0.1e1 / 0.2e1, alpha2) - (beta1 + 4 * c + sigma + 2) / beta1 * besselk(0.3e1 / 0.2e1, alpha2));
a25=(alpha1 * (beta1 + 2) / beta1 * besseli(0.1e1 / 0.2e1, alpha1) - (beta1 + 4 * c + sigma + 2) / beta1 * besseli(0.3e1 / 0.2e1, alpha1)) ;
a26=(alpha2 * (beta1 + 2) / beta1 * besseli(0.1e1 / 0.2e1, alpha2) - (beta1 + 4 * c + sigma + 2) / beta1 * besseli(0.3e1 / 0.2e1, alpha2));
%a31=0;a32=0;
a33= ((alpha1 ^ 2 * beta2 * c * s1 * s2 + alpha1 ^ 2 * beta2 * s1 * s2 - alpha1 ^ 2 * c * s1 - beta2 * s1 * s2 * sigma + alpha1 ^ 2 * c - alpha1 ^ 2 * s1 + alpha1 ^ 2 + s1 * sigma - sigma) / beta2 / s1 / c / s2 * besselk(0.3e1 / 0.2e1, alpha1) / 0.2e1 + 0.1e1 / beta2 / s2 / c * (alpha1 ^ 2 * c + alpha1 ^ 2 - sigma) * besselk(0.5e1 / 0.2e1, alpha1) * alpha1 / 0.2e1) ;
a34= ((alpha2 ^ 2 * beta2 * c * s1 * s2 + alpha2 ^ 2 * beta2 * s1 * s2 - alpha2 ^ 2 * c * s1 - beta2 * s1 * s2 * sigma + alpha2 ^ 2 * c - alpha2 ^ 2 * s1 + alpha2 ^ 2 + s1 * sigma - sigma) / beta2 / s1 / c / s2 * besselk(0.3e1 / 0.2e1, alpha2) / 0.2e1 + 0.1e1 / beta2 / s2 / c * (alpha2 ^ 2 * c + alpha2 ^ 2 - sigma) * besselk(0.5e1 / 0.2e1, alpha2) * alpha2 / 0.2e1) ;
a35= ((alpha1 ^ 2 * beta2 * c * s1 * s2 + alpha1 ^ 2 * beta2 * s1 * s2 - alpha1 ^ 2 * c * s1 - beta2 * s1 * s2 * sigma + alpha1 ^ 2 * c - alpha1 ^ 2 * s1 + alpha1 ^ 2 + s1 * sigma - sigma) / beta2 / s1 / c / s2 * besseli(0.3e1 / 0.2e1, alpha1) / 0.2e1 - 0.1e1 / beta2 / s2 / c * (alpha1 ^ 2 * c + alpha1 ^ 2 - sigma) * alpha1 * besseli(0.5e1 / 0.2e1, alpha1) / 0.2e1) ;
a36= ((alpha2 ^ 2 * beta2 * c * s1 * s2 + alpha2 ^ 2 * beta2 * s1 * s2 - alpha2 ^ 2 * c * s1 - beta2 * s1 * s2 * sigma + alpha2 ^ 2 * c - alpha2 ^ 2 * s1 + alpha2 ^ 2 + s1 * sigma - sigma) / beta2 / s1 / c / s2 * besseli(0.3e1 / 0.2e1, alpha2) / 0.2e1 - 0.1e1 / beta2 / s2 / c * (alpha2 ^ 2 * c + alpha2 ^ 2 - sigma) * besseli(0.5e1 / 0.2e1, alpha2) * alpha2 / 0.2e1) ;
a41=2 / b ^ 3; a42= 2;
a43= 2 * b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha1) ;
a44= 2 * b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha2);
a45= 2 * b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha1) ;
a46= 2 * b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha2);
a51= -1 / b ^ 3 ;a52= 2 ;
a53= (-b ^ (-0.1e1 / 0.2e1) * alpha1 * besselk(0.1e1 / 0.2e1, b * alpha1) - b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha1)) ;
a54= (-b ^ (-0.1e1 / 0.2e1) * alpha2 * besselk(0.1e1 / 0.2e1, b * alpha2) - b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha2)) ;
a55= (b ^ (-0.1e1 / 0.2e1) * alpha1 * besseli(0.1e1 / 0.2e1, b * alpha1) - b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha1)) ;
a56= (b ^ (-0.1e1 / 0.2e1) * alpha2 * besseli(0.1e1 / 0.2e1, b * alpha2) - b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha2)) ;
%a61=0; a62=0;
a63= 0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha1 ^ 2 * b ^ 2 * c + alpha1 ^ 2 * b ^ 2 - b ^ 2 * sigma) * besselk(0.3e1 / 0.2e1, b * alpha1)/ 0.2e1;
a64= 0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha2 ^ 2 * b ^ 2 * c + alpha2 ^ 2 * b ^ 2 - b ^ 2 * sigma) * besselk(0.3e1 / 0.2e1, b * alpha2) / 0.2e1;
a65= 0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha1 ^ 2 * b ^ 2 * c + alpha1 ^ 2 * b ^ 2 - b ^ 2 * sigma) * besseli(0.3e1 / 0.2e1, b * alpha1) / 0.2e1 ;
a66= 0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha2 ^ 2 * b ^ 2 * c + alpha2 ^ 2 * b ^ 2 - b ^ 2 * sigma) * besseli(0.3e1 / 0.2e1, b * alpha2) / 0.2e1;
% Construct the 6x6 matrix A and column vector B
A = [a11, a12, a13, a14, a15, a16;
a21, a22, a23, a24, a25, a26;
0, 0, a33, a34, a35, a36;
a41, a42, a43, a44, a45, a46;
a51, a52, a53, a54, a55, a56;
0, 0, a63, a64, a65, a66];
B = [1; -1; 0; 0; 0; 0]; % Example: solve Ax=B for unit source
x = A \ B;
x1 = x(1);
val =( 3/(4*pi*sigma^2*(rho1 - 1 - 3*x1)));
end
function val = U_function(sigma, p)
% Unpack parameters
c = p.c; b = p.b; s1 = p.s1; s2 = p.s2;
beta1 = p.beta1; beta2 = p.beta2; j = p.j;
rho1 = p.rho1;
% Calculations
term_alpha = (j*s1*sigma + (sigma + 4*c*s1)/(1+c));
root_val = sqrt(term_alpha^2 - (4*s1*sigma*(j*sigma + 4*c)/(1+c)));
alpha1 = sqrt((term_alpha + root_val) / 2);
alpha2 = sqrt((term_alpha - root_val) / 2);
% ... [Insert your existing matrix assignments a11...a66 here] ...
a11=2; a12=2;
a13=2*besselk(1.5, alpha1); a14=2*besselk(1.5, alpha2);
a15=2*besseli(1.5, alpha1); a16=2*besseli(1.5, alpha2);
%b1=U;
a21= -(4*beta1*c + 2*beta1 + 1);
a22= -2*(2*beta1*c - 2*beta1 - 1);
a23=-alpha1*(2*beta1+1)*besselk(0.5,alpha1)-(4*beta1*c+beta1*sigma+2*beta1+1)*besselk(1.5,alpha1);
a24=-alpha2*(2*beta1+1)*besselk(0.5,alpha2)-(4*beta1*c+beta1*sigma+2*beta1+1)*besselk(1.5,alpha2);
a25= alpha1*(2*beta1+1)*besseli(0.5,alpha1)-(4*beta1*c+beta1*sigma+2*beta1+1)*besseli(1.5,alpha1);
a26= alpha2*(2*beta1+1)*besseli(0.5,alpha2)-(4*beta1*c+beta1*sigma+2*beta1+1)*besseli(1.5,alpha2);
%b2=-U;
%a31=0; a32=0;
a33=-(alpha1^2*beta2*c*s2-alpha1^2*c*s1*s2-alpha1^2*beta2*c+alpha1^2*beta2*s2-alpha1^2*s1*s2-alpha1^2*beta2-beta2*s2*sigma+s1*s2*sigma+beta2*sigma)/s2/c/s1*besselk(1.5,alpha1)/2+beta2/s1/c*(alpha1^2*c+alpha1^2-sigma)* besselk(2.5,alpha1)*alpha1/2;
a34=-(alpha2^2*beta2*c*s2-alpha2^2*c*s1*s2-alpha2^2*beta2*c+alpha2^2*beta2*s2-alpha2^2*s1*s2-alpha2^2*beta2-beta2*s2*sigma+s1*s2*sigma+beta2*sigma)/s2/c/s1*besselk(1.5,alpha2)/2+beta2/s1/c*(alpha2^2*c+alpha2^2-sigma)* besselk(2.5,alpha2)*alpha2/2;
a35=-(alpha1^2*beta2*c*s2-alpha1^2*c*s1*s2-alpha1^2*beta2*c+alpha1^2*beta2*s2-alpha1^2*s1*s2-alpha1^2*beta2-beta2*s2*sigma+s1*s2*sigma+beta2*sigma)/s2/c/s1*besseli(1.5,alpha1)/2-beta2/s1/c*(alpha1^2*c+alpha1^2-sigma)*alpha1*besseli(2.5,alpha1) /2;
a36=-(alpha2^2*beta2*c*s2-alpha2^2*c*s1*s2-alpha2^2*beta2*c+alpha2^2*beta2*s2-alpha2^2*s1*s2-alpha2^2*beta2-beta2*s2*sigma+s1*s2*sigma+beta2*sigma)/s2/c/s1*besseli(1.5,alpha2)/2-beta2/s1/c*(alpha2^2*c+alpha2^2-sigma)* besseli(2.5,alpha2)*alpha2/2;
a41=2 / b^3; a42= 2;
a43=2 * b^(-1.5) * besselk(1.5, b * alpha1);
a44=2 * b^(-1.5) * besselk(1.5, b * alpha2);
a45=2 * b^(-1.5) * besseli(1.5, b * alpha1);
a46=2 * b^(-1.5) * besseli(1.5, b * alpha2);
a51=-1/b^3; a52= 2;
a53= (-alpha1 * b^(-0.5) * besselk(0.5, b * alpha1) - b^(-1.5) * besselk(1.5, b * alpha1));
a54= (-alpha2 * b^(-0.5) * besselk(0.5, b * alpha2) - b^(-1.5) * besselk(1.5, b * alpha2));
a55= (alpha1 * b^(-0.5) * besseli(0.5, b * alpha1) - b^(-1.5) * besseli(1.5, b * alpha1));
a56= (alpha2 * b^(-0.5) * besseli(0.5, b * alpha2) - b^(-1.5) * besseli(1.5, b * alpha2));
%a61=0; a62=0;
a63= b^(-0.5) * (alpha1^2 * c + alpha1^2 - sigma) / c * besselk(1.5, b * alpha1) / 0.2e1;
a64= b^(-0.5) * (alpha2^2 * c + alpha2^2 - sigma) / c * besselk(1.5, b * alpha2) / 0.2e1;
a65= b^(-0.5) * (alpha1^2 * c + alpha1^2 - sigma) / c * besseli(1.5, b * alpha1) / 0.2e1;
a66= b^(-0.5) * (alpha2^2 * c + alpha2^2 - sigma) / c * besseli(1.5, b * alpha2) / 0.2e1;
% Numerator/Denominator calculation (truncated for readability)
nm1= a12*a33*a45*a56*a64+a12*a33*a46*a54*a65-a12*a33*a46*a55*a64-a12*a34*a43*a55*a66+a12*a34*a43*a56*a65+a12*a34*a45*a53*a66-a12*a34*a45*a56*a63-a12*a34 * a46 * a53 * a65 + a12 * a34 * a46 * a55 * a63 - a12 * a35 * a43 * a56 * a64 + a12 * a35 * a44 * a56 * a63 + a12 * a35 * a46 * a53 * a64 - a12 * a35 * a46 * a54 * a63 + a12 * a36 * a43 * a55 * a64 - a12 * a36 * a44 * a55 * a63 - a12 * a36 * a45 * a53 * a64 + a12 * a36 * a45 * a54 * a63 + a14 * a35 * a42 * a53 * a66 - a13 * a35 * a42 * a54 * a66 + a13 * a35 * a44 * a52 * a66 + a13 * a36 * a42 * a54 * a65 - a13 * a36 * a44 * a52 * a65 - a14 * a35 * a43 * a52 * a66 - a14 * a36 * a42 * a53 * a65 + a14 * a36 * a43 * a52 * a65 + a12 * a35 * a43 * a54 * a66 - a12 * a35 * a44 * a53 * a66 - a12 * a36 * a43 * a54 * a65 + a12 * a36 * a44 * a53 * a65 - a16 * a34 * a42 * a55 * a63 - a16 * a35 * a42 * a53 * a64 - a16 * a34 * a43 * a52 * a65 + a16 * a34 * a45 * a52 * a63 + a16 * a35 * a42 * a54 * a63 + a16 * a35 * a43 * a52 * a64 - a16 * a35 * a44 * a52 * a63 + a15 * a33 * a42 * a54 * a66 - a15 * a33 * a42 * a56 * a64 - a15 * a33 * a44 * a52 * a66 + a15 * a33 * a46 * a52 * a64 - a15 * a34 * a42 * a53 * a66 + a15 * a34 * a42 * a56 * a63 + a15 * a34 * a43 * a52 * a66 - a15 * a34 * a46 * a52 * a63 + a15 * a36 * a42 * a53 * a64 - a15 * a36 * a42 * a54 * a63 - a15 * a36 * a43 * a52 * a64 + a15 * a36 * a44 * a52 * a63 - a16 * a33 * a42 * a54 * a65 - a14 * a33 * a42 * a55 * a66 + a14 * a33 * a42 * a56 * a65 + a14 * a33 * a45 * a52 * a66 - a14 * a33 * a46 * a52 * a65 - a14 * a35 * a42 * a56 * a63 + a14 * a35 * a46 * a52 * a63 + a14 * a36 * a42 * a55 * a63 - a14 * a36 * a45 * a52 * a63 + a13 * a34 * a42 * a55 * a66 - a13 * a34 * a42 * a56 * a65 - a13 * a34 * a45 * a52 * a66 + a13 * a34 * a46 * a52 * a65 + a13 * a35 * a42 * a56 * a64 - a13 * a35 * a46 * a52 * a64 - a13 * a36 * a42 * a55 * a64 + a13 * a36 * a45 * a52 * a64 + a12 * a33 * a44 * a55 * a66 - a12 * a33 * a44 * a56 * a65 - a12 * a33 * a45 * a54 * a66 + a16 * a33 * a42 * a55 * a64 + a16 * a33 * a44 * a52 * a65 - a16 * a33 * a45 * a52 * a64 + a16 * a34 * a42 * a53 * a65 - a25 * a34 * a46 * a52 * a63 + a25 * a36 * a42 * a53 * a64 - a25 * a36 * a42 * a54 * a63 + a22 * a33 * a44 * a55 * a66 - a22 * a33 * a44 * a56 * a65 - a22 * a33 * a45 * a54 * a66 + a22 * a33 * a45 * a56 * a64 + a22 * a33 * a46 * a54 * a65 - a22 * a33 * a46 * a55 * a64 - a22 * a34 * a43 * a55 * a66 + a22 * a34 * a43 * a56 * a65 + a22 * a34 * a45 * a53 * a66 - a22 * a34 * a45 * a56 * a63 - a22 * a34 * a46 * a53 * a65 + a22 * a34 * a46 * a55 * a63 - a22 * a35 * a43 * a56 * a64 + a22 * a35 * a44 * a56 * a63 + a22 * a35 * a46 * a53 * a64 - a22 * a35 * a46 * a54 * a63 - a22 * a36 * a43 * a54 * a65 + a22 * a36 * a44 * a53 * a65 + a22 * a35 * a43 * a54 * a66 - a23 * a36 * a44 * a52 * a65 - a23 * a35 * a42 * a54 * a66 + a23 * a35 * a44 * a52 * a66 + a23 * a36 * a42 * a54 * a65 + a24 * a35 * a42 * a53 * a66 - a24 * a35 * a43 * a52 * a66 - a24 * a36 * a42 * a53 * a65 + a24 * a36 * a43 * a52 * a65 - a22 * a35 * a44 * a53 * a66 - a26 * a33 * a42 * a54 * a65 + a26 * a33 * a42 * a55 * a64 + a26 * a33 * a44 * a52 * a65 - a26 * a34 * a42 * a55 * a63 - a26 * a33 * a45 * a52 * a64 + a26 * a34 * a42 * a53 * a65 - a26 * a34 * a43 * a52 * a65 + a26 * a34 * a45 * a52 * a63 - a26 * a35 * a42 * a53 * a64 + a26 * a35 * a42 * a54 * a63 + a26 * a35 * a43 * a52 * a64 - a26 * a35 * a44 * a52 * a63 + a22 * a36 * a43 * a55 * a64 - a22 * a36 * a44 * a55 * a63 - a22 * a36 * a45 * a53 * a64 + a22 * a36 * a45 * a54 * a63 + a23 * a34 * a42 * a55 * a66 - a23 * a34 * a42 * a56 * a65 - a23 * a34 * a45 * a52 * a66 + a23 * a34 * a46 * a52 * a65 + a23 * a35 * a42 * a56 * a64 - a23 * a35 * a46 * a52 * a64 - a23 * a36 * a42 * a55 * a64 + a23 * a36 * a45 * a52 * a64 - a24 * a33 * a42 * a55 * a66 + a24 * a33 * a42 * a56 * a65 + a24 * a33 * a45 * a52 * a66 - a24 * a33 * a46 * a52 * a65 - a24 * a35 * a42 * a56 * a63 + a24 * a35 * a46 * a52 * a63 + a24 * a36 * a42 * a55 * a63 - a24 * a36 * a45 * a52 * a63 + a25 * a33 * a42 * a54 * a66 - a25 * a33 * a42 * a56 * a64 - a25 * a33 * a44 * a52 * a66 + a25 * a33 * a46 * a52 * a64 - a25 * a34 * a42 * a53 * a66 + a25 * a34 * a42 * a56 * a63 + a25 * a34 * a43 * a52 * a66 - a25 * a36 * a43 * a52 * a64 + a25 * a36 * a44 * a52 * a63;
nm2=-a16*a21*a33*a42*a55*a64-a16*a21*a33*a44*a52*a65+a16*a21*a33*a45*a52*a64-a16*a21*a34*a42*a53*a65+a16*a21*a34*a42*a55*a63+a16*a21*a34*a43*a52*a65-a16 * a21 * a34 * a45 * a52 * a63 + a16 * a21 * a35 * a42 * a53 * a64 - a16 * a21 * a35 * a42 * a54 * a63 - a16 * a21 * a35 * a43 * a52 * a64 + a16 * a21 * a35 * a44 * a52 * a63 - a16 * a22 * a33 * a41 * a54 * a65 + a16 * a22 * a33 * a41 * a55 * a64 + a16 * a22 * a33 * a44 * a51 * a65 - a16 * a22 * a33 * a45 * a51 * a64 + a16 * a22 * a34 * a41 * a53 * a65 - a16 * a22 * a34 * a41 * a55 * a63 - a16 * a22 * a34 * a43 * a51 * a65 + a16 * a22 * a34 * a45 * a51 * a63 - a16 * a22 * a35 * a41 * a53 * a64 + a16 * a22 * a35 * a41 * a54 * a63 + a16 * a22 * a35 * a43 * a51 * a64 - a16 * a22 * a35 * a44 * a51 * a63 - a16 * a23 * a34 * a41 * a52 * a65 + a16 * a23 * a34 * a42 * a51 * a65 + a16 * a23 * a35 * a41 * a52 * a64 - a16 * a23 * a35 * a42 * a51 * a64 + a16 * a24 * a33 * a41 * a52 * a65 - a16 * a24 * a33 * a42 * a51 * a65 - a16 * a24 * a35 * a41 * a52 * a63 + a16 * a24 * a35 * a42 * a51 * a63 - a16 * a25 * a33 * a41 * a52 * a64 + a16 * a25 * a33 * a42 * a51 * a64 + a16 * a25 * a34 * a41 * a52 * a63 - a16 * a25 * a34 * a42 * a51 * a63 + a14 * a26 * a35 * a41 * a52 * a63 - a14 * a26 * a35 * a42 * a51 * a63 - a15 * a21 * a33 * a42 * a54 * a66 + a15 * a21 * a33 * a42 * a56 * a64 + a15 * a21 * a33 * a44 * a52 * a66 - a15 * a21 * a33 * a46 * a52 * a64 + a15 * a21 * a34 * a42 * a53 * a66 - a15 * a21 * a34 * a42 * a56 * a63 - a15 * a21 * a34 * a43 * a52 * a66 + a15 * a21 * a34 * a46 * a52 * a63 - a15 * a21 * a36 * a42 * a53 * a64 + a15 * a21 * a36 * a42 * a54 * a63 + a15 * a21 * a36 * a43 * a52 * a64 - a15 * a21 * a36 * a44 * a52 * a63 + a15 * a22 * a33 * a41 * a54 * a66 - a15 * a22 * a33 * a41 * a56 * a64 - a15 * a22 * a33 * a44 * a51 * a66 + a15 * a22 * a33 * a46 * a51 * a64 - a15 * a22 * a34 * a41 * a53 * a66 + a15 * a22 * a34 * a41 * a56 * a63 + a15 * a22 * a34 * a43 * a51 * a66 - a15 * a22 * a34 * a46 * a51 * a63 + a15 * a22 * a36 * a41 * a53 * a64 - a15 * a22 * a36 * a41 * a54 * a63 - a15 * a22 * a36 * a43 * a51 * a64 + a15 * a22 * a36 * a44 * a51 * a63 + a15 * a23 * a34 * a41 * a52 * a66 - a15 * a23 * a34 * a42 * a51 * a66 - a15 * a23 * a36 * a41 * a52 * a64 + a15 * a23 * a36 * a42 * a51 * a64 - a15 * a24 * a33 * a41 * a52 * a66 + a15 * a24 * a33 * a42 * a51 * a66 + a15 * a24 * a36 * a41 * a52 * a63 - a15 * a24 * a36 * a42 * a51 * a63 + a15 * a26 * a33 * a41 * a52 * a64 - a15 * a26 * a33 * a42 * a51 * a64 - a15 * a26 * a34 * a41 * a52 * a63 + a15 * a26 * a34 * a42 * a51 * a63 + a16 * a21 * a33 * a42 * a54 * a65 - a13 * a25 * a36 * a42 * a51 * a64 + a13 * a26 * a34 * a41 * a52 * a65 - a13 * a26 * a34 * a42 * a51 * a65 - a13 * a26 * a35 * a41 * a52 * a64 + a13 * a26 * a35 * a42 * a51 * a64 + a14 * a21 * a33 * a42 * a55 * a66 - a14 * a21 * a33 * a42 * a56 * a65 - a14 * a21 * a33 * a45 * a52 * a66 + a14 * a21 * a33 * a46 * a52 * a65 + a14 * a21 * a35 * a42 * a56 * a63 - a14 * a21 * a35 * a46 * a52 * a63 - a14 * a21 * a36 * a42 * a55 * a63 + a14 * a21 * a36 * a45 * a52 * a63 - a14 * a22 * a33 * a41 * a55 * a66 + a14 * a22 * a33 * a41 * a56 * a65 + a14 * a22 * a33 * a45 * a51 * a66 - a14 * a22 * a33 * a46 * a51 * a65 - a14 * a22 * a35 * a41 * a56 * a63 + a14 * a22 * a35 * a46 * a51 * a63 + a14 * a22 * a36 * a41 * a55 * a63 - a14 * a22 * a36 * a45 * a51 * a63 + a14 * a25 * a33 * a41 * a52 * a66 - a14 * a25 * a33 * a42 * a51 * a66 - a14 * a25 * a36 * a41 * a52 * a63 + a14 * a25 * a36 * a42 * a51 * a63 - a14 * a26 * a33 * a41 * a52 * a65 + a14 * a26 * a33 * a42 * a51 * a65 - a12 * a26 * a34 * a41 * a53 * a65 + a12 * a26 * a34 * a41 * a55 * a63 + a12 * a26 * a34 * a43 * a51 * a65 - a12 * a26 * a34 * a45 * a51 * a63 + a12 * a26 * a35 * a41 * a53 * a64 - a12 * a26 * a35 * a41 * a54 * a63 - a12 * a26 * a35 * a43 * a51 * a64 + a12 * a26 * a35 * a44 * a51 * a63 - a13 * a21 * a34 * a42 * a55 * a66 + a13 * a21 * a34 * a42 * a56 * a65 + a13 * a21 * a34 * a45 * a52 * a66 - a13 * a21 * a34 * a46 * a52 * a65 - a13 * a21 * a35 * a42 * a56 * a64 + a13 * a21 * a35 * a46 * a52 * a64 + a13 * a21 * a36 * a42 * a55 * a64 - a13 * a21 * a36 * a45 * a52 * a64 + a13 * a22 * a34 * a41 * a55 * a66 - a13 * a22 * a34 * a41 * a56 * a65 - a13 * a22 * a34 * a45 * a51 * a66 + a13 * a22 * a34 * a46 * a51 * a65 + a13 * a22 * a35 * a41 * a56 * a64 - a13 * a22 * a35 * a46 * a51 * a64 - a13 * a22 * a36 * a41 * a55 * a64 + a13 * a22 * a36 * a45 * a51 * a64 - a13 * a25 * a34 * a41 * a52 * a66 + a13 * a25 * a34 * a42 * a51 * a66 + a13 * a25 * a36 * a41 * a52 * a64 + a12 * a23 * a34 * a41 * a56 * a65 + a12 * a23 * a34 * a45 * a51 * a66 - a12 * a23 * a34 * a46 * a51 * a65 - a12 * a23 * a35 * a41 * a56 * a64 + a12 * a23 * a35 * a46 * a51 * a64 + a12 * a23 * a36 * a41 * a55 * a64 - a12 * a23 * a36 * a45 * a51 * a64 + a12 * a24 * a33 * a41 * a55 * a66 - a12 * a24 * a33 * a41 * a56 * a65 - a12 * a24 * a33 * a45 * a51 * a66 + a12 * a24 * a33 * a46 * a51 * a65 + a12 * a24 * a35 * a41 * a56 * a63 - a12 * a24 * a35 * a46 * a51 * a63 - a12 * a24 * a36 * a41 * a55 * a63 + a12 * a24 * a36 * a45 * a51 * a63 - a12 * a25 * a33 * a41 * a54 * a66 + a12 * a25 * a33 * a41 * a56 * a64 + a12 * a25 * a33 * a44 * a51 * a66 - a12 * a25 * a33 * a46 * a51 * a64 + a12 * a25 * a34 * a41 * a53 * a66 - a12 * a25 * a34 * a41 * a56 * a63 - a12 * a25 * a34 * a43 * a51 * a66 + a12 * a25 * a34 * a46 * a51 * a63 - a12 * a25 * a36 * a41 * a53 * a64 + a12 * a25 * a36 * a41 * a54 * a63 + a12 * a25 * a36 * a43 * a51 * a64 - a12 * a25 * a36 * a44 * a51 * a63 + a12 * a26 * a33 * a41 * a54 * a65 - a12 * a26 * a33 * a41 * a55 * a64 - a12 * a26 * a33 * a44 * a51 * a65 + a12 * a26 * a33 * a45 * a51 * a64 - a11 * a25 * a36 * a43 * a52 * a64 + a11 * a25 * a36 * a44 * a52 * a63 - a11 * a26 * a33 * a42 * a54 * a65 + a11 * a26 * a33 * a42 * a55 * a64 + a11 * a26 * a33 * a44 * a52 * a65 - a11 * a26 * a33 * a45 * a52 * a64 + a11 * a26 * a34 * a42 * a53 * a65 - a11 * a26 * a34 * a42 * a55 * a63 - a11 * a26 * a34 * a43 * a52 * a65 + a11 * a26 * a34 * a45 * a52 * a63 - a11 * a26 * a35 * a42 * a53 * a64 + a11 * a26 * a35 * a42 * a54 * a63 + a11 * a26 * a35 * a43 * a52 * a64 - a11 * a26 * a35 * a44 * a52 * a63 - a12 * a21 * a33 * a44 * a55 * a66 + a12 * a21 * a33 * a44 * a56 * a65 + a12 * a21 * a33 * a45 * a54 * a66 - a12 * a21 * a33 * a45 * a56 * a64 - a12 * a21 * a33 * a46 * a54 * a65 + a12 * a21 * a33 * a46 * a55 * a64 + a12 * a21 * a34 * a43 * a55 * a66 - a12 * a21 * a34 * a43 * a56 * a65 - a12 * a21 * a34 * a45 * a53 * a66 + a12 * a21 * a34 * a45 * a56 * a63 + a12 * a21 * a34 * a46 * a53 * a65 - a12 * a21 * a34 * a46 * a55 * a63 + a12 * a21 * a35 * a43 * a56 * a64 - a12 * a21 * a35 * a44 * a56 * a63 - a12 * a21 * a35 * a46 * a53 * a64 + a12 * a21 * a35 * a46 * a54 * a63 - a12 * a21 * a36 * a43 * a55 * a64 + a12 * a21 * a36 * a44 * a55 * a63 + a12 * a21 * a36 * a45 * a53 * a64 - a12 * a21 * a36 * a45 * a54 * a63 - a12 * a23 * a34 * a41 * a55 * a66 + a11 * a22 * a36 * a43 * a55 * a64 - a11 * a22 * a36 * a44 * a55 * a63 - a11 * a22 * a36 * a45 * a53 * a64 + a11 * a22 * a36 * a45 * a54 * a63 + a11 * a23 * a34 * a42 * a55 * a66 - a11 * a23 * a34 * a42 * a56 * a65 - a11 * a23 * a34 * a45 * a52 * a66 + a11 * a23 * a34 * a46 * a52 * a65 + a11 * a23 * a35 * a42 * a56 * a64 - a11 * a23 * a35 * a46 * a52 * a64 - a11 * a23 * a36 * a42 * a55 * a64 + a11 * a23 * a36 * a45 * a52 * a64 - a11 * a24 * a33 * a42 * a55 * a66 + a11 * a24 * a33 * a42 * a56 * a65 + a11 * a24 * a33 * a45 * a52 * a66 - a11 * a24 * a33 * a46 * a52 * a65 - a11 * a24 * a35 * a42 * a56 * a63 + a11 * a24 * a35 * a46 * a52 * a63 + a11 * a24 * a36 * a42 * a55 * a63 - a11 * a24 * a36 * a45 * a52 * a63 + a11 * a25 * a33 * a42 * a54 * a66 - a11 * a25 * a33 * a42 * a56 * a64 - a11 * a25 * a33 * a44 * a52 * a66 + a11 * a25 * a33 * a46 * a52 * a64 - a11 * a25 * a34 * a42 * a53 * a66 + a11 * a25 * a34 * a42 * a56 * a63 + a11 * a25 * a34 * a43 * a52 * a66 - a11 * a25 * a34 * a46 * a52 * a63 + a11 * a25 * a36 * a42 * a53 * a64 - a11 * a25 * a36 * a42 * a54 * a63 - a14 * a23 * a36 * a42 * a51 * a65 - a14 * a23 * a35 * a41 * a52 * a66 + a14 * a23 * a35 * a42 * a51 * a66 + a14 * a23 * a36 * a41 * a52 * a65 - a14 * a21 * a35 * a42 * a53 * a66 + a11 * a22 * a33 * a44 * a55 * a66 - a11 * a22 * a33 * a44 * a56 * a65 - a11 * a22 * a33 * a45 * a54 * a66 + a11 * a22 * a33 * a45 * a56 * a64 + a11 * a22 * a33 * a46 * a54 * a65 - a11 * a22 * a33 * a46 * a55 * a64 - a11 * a22 * a34 * a43 * a55 * a66 + a11 * a22 * a34 * a43 * a56 * a65 + a11 * a22 * a34 * a45 * a53 * a66 - a11 * a22 * a34 * a45 * a56 * a63 - a11 * a22 * a34 * a46 * a53 * a65 + a11 * a22 * a34 * a46 * a55 * a63 - a11 * a22 * a35 * a43 * a56 * a64 + a11 * a22 * a35 * a44 * a56 * a63 + a11 * a22 * a35 * a46 * a53 * a64 - a11 * a22 * a35 * a46 * a54 * a63 - a12 * a24 * a35 * a41 * a53 * a66 + a12 * a24 * a35 * a43 * a51 * a66 + a12 * a24 * a36 * a41 * a53 * a65 - a13 * a22 * a35 * a41 * a54 * a66 + a13 * a22 * a35 * a44 * a51 * a66 + a13 * a22 * a36 * a41 * a54 * a65 - a13 * a22 * a36 * a44 * a51 * a65 - a12 * a24 * a36 * a43 * a51 * a65 + a13 * a24 * a36 * a42 * a51 * a65 + a13 * a21 * a35 * a42 * a54 * a66 - a13 * a21 * a35 * a44 * a52 * a66 - a13 * a21 * a36 * a42 * a54 * a65 + a13 * a21 * a36 * a44 * a52 * a65 + a14 * a21 * a35 * a43 * a52 * a66 + a14 * a21 * a36 * a42 * a53 * a65 - a14 * a21 * a36 * a43 * a52 * a65 + a14 * a22 * a35 * a41 * a53 * a66 - a14 * a22 * a35 * a43 * a51 * a66 - a14 * a22 * a36 * a41 * a53 * a65 + a14 * a22 * a36 * a43 * a51 * a65 + a13 * a24 * a35 * a41 * a52 * a66 - a13 * a24 * a35 * a42 * a51 * a66 - a13 * a24 * a36 * a41 * a52 * a65 - a11 * a22 * a36 * a43 * a54 * a65 + a11 * a22 * a36 * a44 * a53 * a65 + a11 * a22 * a35 * a43 * a54 * a66 - a11 * a23 * a36 * a44 * a52 * a65 - a11 * a23 * a35 * a42 * a54 * a66 + a11 * a23 * a35 * a44 * a52 * a66 + a11 * a23 * a36 * a42 * a54 * a65 + a11 * a24 * a35 * a42 * a53 * a66 - a11 * a24 * a35 * a43 * a52 * a66 - a11 * a24 * a36 * a42 * a53 * a65 + a11 * a24 * a36 * a43 * a52 * a65 - a12 * a21 * a35 * a43 * a54 * a66 + a12 * a21 * a35 * a44 * a53 * a66 + a12 * a21 * a36 * a43 * a54 * a65 - a12 * a21 * a36 * a44 * a53 * a65 + a12 * a23 * a35 * a41 * a54 * a66 - a12 * a23 * a35 * a44 * a51 * a66 - a12 * a23 * a36 * a41 * a54 * a65 + a12 * a23 * a36 * a44 * a51 * a65 - a11 * a22 * a35 * a44 * a53 * a66;
A1=nm1/nm2; %without u(sigma)
l = sqrt(4 * c * s1 / (1 + c));j1 = 2 + (s1/beta2) + 2*(s1 / s2);
FA=(6*pi*(c+1)*(2*beta1+1)*(l^2+l*j1+j1))/((3*beta1+1)*(c+1)*(l^2+l*j1+j1)-j1*c*(2*beta1+1));
val=3*FA/((4*pi*sigma^2)*(rho1-1-3*A1));
end
function f = talbot_inversion(F, t)
N = 25;
f = 0;
% Talbot Contour
for k = 1:N
theta = (k-0.5)*pi/N;
s = (N/t) * (theta*cot(theta) + 1i*theta);
ds = (N/t) * (theta*(-csc(theta)^2) + cot(theta) + 1i);
% Talbot weight
f = f + exp(s*t) * F(s) * ds;
end
f = real(f / (2*pi*1i));
end
Torsten
on 8 Jul 2026 at 17:17
If @Shreen El-Sapa simply tries to compute lim(t->Inf) u_s1(t), it's of course possible to compute this limit the way you do. But I wonder why an extra function U_func_steady is needed for this and how this function is obtained.
Shreen El-Sapa
on 8 Jul 2026 at 17:56
Shreen El-Sapa
on 8 Jul 2026 at 17:57
Edited: Torsten
on 8 Jul 2026 at 18:54
Try this, but the limit case does not appear. At t tends to infity the curves should go to 1:
warning('off')
test_1()
function test_1()
clc; close all;
t_values = logspace(0, 2, 100);
rho1_vals = [0, 10, 50, 100]; % 4 specific curves to compare
params = struct('s1',1, 's2',100, 'c',0.000001,'beta1',1000000, 'beta2',1000000, 'b',2,'j',0.1);
figure; hold on;
colors = ['r', 'g', 'b', 'm']; % Three colors for the three curves
for k = 1:length(rho1_vals)
params.rho1 = rho1_vals(k);
u_s1 = zeros(size(t_values));
u_steady = zeros(size(t_values));
for i = 1:length(t_values)
% Compute both inversions at the current beta1
u_s1(i) = talbot_inversion(@(s) U_func_s1(s, params), t_values(i));
u_steady(i) = talbot_inversion(@(s) U_func_steady(params), t_values(i));
end
% Calculate ratio
ratio = u_s1 ./ (u_steady + eps);
% Plot
plot(t_values, ratio, 'Color', colors(k), 'LineWidth', 2, ...
'DisplayName', sprintf('\\rho=%.2f', rho1_vals(k)));
end
set(gca, 'XScale', 'log');
xlabel('Time (t)'); ylabel('Ratio (U_{uns}/U_{std})');
legend('Location', 'best'); grid on;
% title('Effect of Micropolar Coupling on Dispersion Ratio');
end
% --- Rename your existing logic to avoid conflicts ---
function val = U_func_s1(sigma, p)
% Copy your logic from the original calc_s1.m U_function here
% Ensure it returns 'val'
% Unpack
b = p.b; c = p.c; s1 = p.s1; s2 = p.s2;
beta1 = p.beta1; beta2 = p.beta2; j = p.j; rho1 = p.rho1;
% Define S (sum) and P (product)
S = (4*c*s1 + (1 + j*s1*(c + 1))*sigma) / (c + 1);
P = (sigma*s1*(4*c + j*sigma)) / (c + 1);
% Calculate alpha1 and alpha2
alpha1 = sqrt((S - sqrt(S^2 - 4*P)) / 2);
alpha2 = sqrt((S + sqrt(S^2 - 4*P)) / 2);
%A1,B1,C1,D1,E1,F1
a11= -2 ;a12= - 2 ;
a13=- 2 *besselk(0.3e1 / 0.2e1, alpha1) ;
a14=- 2 *besselk(0.3e1 / 0.2e1, alpha2) ;
a15=- 2 *besseli(0.3e1 / 0.2e1, alpha1) ;
a16=- 2 *besseli(0.3e1 / 0.2e1, alpha2);
a21= -(beta1 + 6) / beta1 ;
a22= 2;
a23= (-alpha1 * (beta1 + 2) / beta1 * besselk(0.1e1 / 0.2e1, alpha1) - (beta1 + sigma + 6) / beta1 * besselk(0.3e1 / 0.2e1, alpha1)) ;
a24= (-alpha2 * (beta1 + 2) / beta1 * besselk(0.1e1 / 0.2e1, alpha2) - (beta1 + sigma + 6) / beta1 * besselk(0.3e1 / 0.2e1, alpha2)) ;
a25= (alpha1 * (beta1 + 2) / beta1 * besseli(0.1e1 / 0.2e1, alpha1) - (beta1 + sigma + 6) / beta1 * besseli(0.3e1 / 0.2e1, alpha1)) ;
a26= (alpha2 * (beta1 + 2) / beta1 * besseli(0.1e1 / 0.2e1, alpha2) - (beta1 + sigma + 6) / beta1 * besseli(0.3e1 / 0.2e1, alpha2)) ;
%a31=0;a32=0;
a33=((alpha1 ^ 2 * beta2 * c * s1 * s2 + alpha1 ^ 2 * beta2 * s1 * s2 - alpha1 ^ 2 * c * s1 - beta2 * s1 * s2 * sigma + alpha1 ^ 2 * c - alpha1 ^ 2 * s1 + alpha1 ^ 2 + s1 * sigma - sigma) / beta2 / s1 / c / s2 * besselk(0.3e1 / 0.2e1, alpha1) / 0.2e1 + 0.1e1 / beta2 / s2 / c * (alpha1 ^ 2 * c + alpha1 ^ 2 - sigma) * besselk(0.5e1 / 0.2e1, alpha1) * alpha1 / 0.2e1) ;
a34=((alpha2 ^ 2 * beta2 * c * s1 * s2 + alpha2 ^ 2 * beta2 * s1 * s2 - alpha2 ^ 2 * c * s1 - beta2 * s1 * s2 * sigma + alpha2 ^ 2 * c - alpha2 ^ 2 * s1 + alpha2 ^ 2 + s1 * sigma - sigma) / beta2 / s1 / c / s2 * besselk(0.3e1 / 0.2e1, alpha2) / 0.2e1 + 0.1e1 / beta2 / s2 / c * (alpha2 ^ 2 * c + alpha2 ^ 2 - sigma) * besselk(0.5e1 / 0.2e1, alpha2) * alpha2 / 0.2e1) ;
a35=((alpha1 ^ 2 * beta2 * c * s1 * s2 + alpha1 ^ 2 * beta2 * s1 * s2 - alpha1 ^ 2 * c * s1 - beta2 * s1 * s2 * sigma + alpha1 ^ 2 * c - alpha1 ^ 2 * s1 + alpha1 ^ 2 + s1 * sigma - sigma) / beta2 / s1 / c / s2 * besseli(0.3e1 / 0.2e1, alpha1) / 0.2e1 - 0.1e1 / beta2 / s2 / c * (alpha1 ^ 2 * c + alpha1 ^ 2 - sigma) * alpha1 * besseli(0.5e1 / 0.2e1, alpha1) / 0.2e1) ;
a36=((alpha2 ^ 2 * beta2 * c * s1 * s2 + alpha2 ^ 2 * beta2 * s1 * s2 - alpha2 ^ 2 * c * s1 - beta2 * s1 * s2 * sigma + alpha2 ^ 2 * c - alpha2 ^ 2 * s1 + alpha2 ^ 2 + s1 * sigma - sigma) / beta2 / s1 / c / s2 * besseli(0.3e1 / 0.2e1, alpha2) / 0.2e1 - 0.1e1 / beta2 / s2 / c * (alpha2 ^ 2 * c + alpha2 ^ 2 - sigma) * besseli(0.5e1 / 0.2e1, alpha2) * alpha2 / 0.2e1) ;
a41= -2 / b ^ 3 ;
a42=- 2 ;
a43=- 2 * b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha1) ;
a44=- 2 * b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha2) ;
a45=- 2 * b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha1) ;
a46=- 2 * b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha2);
a51= -1 / b ^ 3 ;
a52= 2 ;
a53= (-b ^ (-0.1e1 / 0.2e1) * alpha1 * besselk(0.1e1 / 0.2e1, b * alpha1) - b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha1)) ;
a54= (-b ^ (-0.1e1 / 0.2e1) * alpha2 * besselk(0.1e1 / 0.2e1, b * alpha2) - b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha2)) ;
a55= (b ^ (-0.1e1 / 0.2e1) * alpha1 * besseli(0.1e1 / 0.2e1, b * alpha1) - b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha1)) ;
a56= (b ^ (-0.1e1 / 0.2e1) * alpha2 * besseli(0.1e1 / 0.2e1, b * alpha2) - b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha2)) ;
%a61 0;a62=0;
a63=0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha1 ^ 2 * b ^ 2 * c + alpha1 ^ 2 * b ^ 2 - b ^ 2 * sigma) * besselk(0.3e1 / 0.2e1, b * alpha1) / 0.2e1 ;
a64=0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha2 ^ 2 * b ^ 2 * c + alpha2 ^ 2 * b ^ 2 - b ^ 2 * sigma) * besselk(0.3e1 / 0.2e1, b * alpha2) / 0.2e1 ;
a65=0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha1 ^ 2 * b ^ 2 * c + alpha1 ^ 2 * b ^ 2 - b ^ 2 * sigma) * besseli(0.3e1 / 0.2e1, b * alpha1)/ 0.2e1 ;
a66=0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha2 ^ 2 * b ^ 2 * c + alpha2 ^ 2 * b ^ 2 - b ^ 2 * sigma) * besseli(0.3e1 / 0.2e1, b * alpha2) / 0.2e1;
% Construct the 6x6 matrix A and column vector B
AA = [a11, a12, a13, a14, a15, a16;
a21, a22, a23, a24, a25, a26;
0, 0, a33, a34, a35, a36;
a41, a42, a43, a44, a45, a46;
a51, a52, a53, a54, a55, a56;
0, 0, a63, a64, a65, a66];
BB = [1; -1; 0; 0; 0; 0]; % Example: solve Ax=B for unit source
x = AA \ BB;
x1 = x(1);
val =( 3/(4*pi*sigma^2*(rho1 - 1 - 3*x1)));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function val0 = U_func_steady(p)
% Copy your logic from the original calc_steady.m U_function here
% Ensure it returns 'val0'
% Unpack
b = p.b; c = p.c; s1 = p.s1; s2 = p.s2;
beta1 = p.beta1; beta2 = p.beta2;
% Define S (sum) and P (product)
alpha = sqrt(4*c*s1 /(c + 1));
a11 = -2 ; a12= - 2 ;a13= - 2 ;a14= - 2 ;
a15=- 2 * besselk(0.3e1 / 0.2e1, alpha) ;
a16=- 2 * besseli(0.3e1 / 0.2e1, alpha) ;
a21 = -(beta1 + 6) / beta1 ;
a22=2 ;
a23=2 * (2 * beta1 - 5 * c - 3) / beta1 ;
a24=(beta1 + 2 * c) / beta1 ;
a25=(-alpha * (beta1 + 2) / beta1 * besselk(0.1e1 / 0.2e1, alpha) - (beta1 + 6) / beta1 * besselk(0.3e1 / 0.2e1, alpha)) ;
a26=(alpha * (beta1 + 2) / beta1 * besseli(0.1e1 / 0.2e1, alpha) - (beta1 + 6) / beta1 * besseli(0.3e1 / 0.2e1, alpha)) ;
%a31=0;a32=0;a33=0;a34=0;
a35= ((1 / beta2 * (1 + c) * alpha ^ 3 / s2 / c * besselk(0.5e1 / 0.2e1, alpha)) / 0.2e1 + ((1 + c) * (alpha ^ 2 * beta2 * s1 * s2 - alpha ^ 2 * s1 + alpha ^ 2) / beta2 / s2 / c / s1 * besselk(0.3e1 / 0.2e1, alpha)) / 0.2e1) ;
a36= (-(1 / beta2 * (1 + c) * alpha ^ 3 / s2 / c * besseli(0.5e1 / 0.2e1, alpha)) / 0.2e1 + ((1 + c) * (alpha ^ 2 * beta2 * s1 * s2 - alpha ^ 2 * s1 + alpha ^ 2) / beta2 / s2 / c / s1 * besseli(0.3e1 / 0.2e1, alpha)) / 0.2e1) ;
a41=-2 / b ^ 3 ;
a42=- 2;
a43=- 2 * b ^ 2 ;
a44=- 2 / b ;
a45=- 2 * b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha) ;
a46=- 2 * b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha);
a51= -1 / b ^ 3 ;
a52= 2 ;
a53=4 * b ^ 2 ;
a54= 1 / b ;
a55=(-b ^ (-0.1e1 / 0.2e1) * alpha * besselk(0.1e1 / 0.2e1, b * alpha) - b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha)) ;
a56=(b ^ (-0.1e1 / 0.2e1) * alpha * besseli(0.1e1 / 0.2e1, b * alpha) - b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha)) ;
%a61=0;a62=0;a63=0;a64=0;
a65= ((1 + c) / c * b ^ (-0.1e1 / 0.2e1) * alpha ^ 2 * besselk(0.3e1 / 0.2e1, b * alpha) ) / 0.2e1;
a66= ((1 + c) / c * b ^ (-0.1e1 / 0.2e1) * alpha ^ 2 * besseli(0.3e1 / 0.2e1, b * alpha)) / 0.2e1;
% Construct the 6x6 matrix A and column vector B
A0 = [a11, a12, a13, a14, a15, a16;
a21, a22, a23, a24, a25, a26;
0, 0, 0, 0, a35, a36;
a41, a42, a43, a44, a45, a46;
a51, a52, a53, a54, a55, a56;
0, 0, 0, 0, a65, a66];
B0 = [1; -1; 0; 0; 0; 0]; % Example: solve Ax=B for unit source
x0 = A0 \ B0;
x0 = x0(1);
% val0 =(3/(4*pi*sigma^2*(1 + 3*x0)));
% val0 =(3/(4*pi*(1 + 3*x0)));
val0 =1/(6*pi*x0);
end
function f = talbot_inversion(F, t)
N = 25;
f = 0;
% Talbot Contour
for k = 1:N
theta = (k-0.5)*pi/N;
s = (N/t) * (theta*cot(theta) + 1i*theta);
ds = (N/t) * (theta*(-csc(theta)^2) + cot(theta) + 1i);
% Talbot weight
f = f + exp(s*t) * F(s) * ds;
end
f = real(f / (2*pi*1i));
end
Torsten
on 8 Jul 2026 at 18:51
If U_steady doesn't depend on t, why do you compute it as
u_steady(i) = talbot_inversion(@(s) U_func_steady(params), t_values(i));
thus as a function of t_values ?
Or do you only need one value, namely
u_steady = talbot_inversion(@(s) U_func_steady(params), Inf);
to compute "ratio" ?
dpb
on 8 Jul 2026 at 18:53
Again, PLEASE FORMAT THE CODE!!!!
_"U_steady does not depend on t"_
Have you verified that? By the way the code is constructed, it is recalculated every timestep and the various terms are, in part, dependent upon the input t.
While it's too dense code to check for sure without more time and patience have; it appears the second function is identical to the first with the exception of some machinations for a val0 variable at the end--no klew is given as to from whence that came, but if it is known what the transfer function should be analytically, then that could simply be coded as it kinda' looks like something similar is being done there in patching in constants.
But, as shown in the above example, the returned value of u_s1(t) doesn't actually have a discernible steady-state value anyway and if there is some analytic form known for which the steady state value is known, then unless the functional can compute that value, something is wrong with it (or the SS value isn't right), but that it doesn't return a plateau is in contradiction to the assumption that it does.
Again, the condition value warnings are probably not a good sign...
Shreen El-Sapa
on 8 Jul 2026 at 19:43
Well, again, you need to look at what you're generating -- plot the two pieces-parts...
The top is the u_s, the middle the SS -- as can be seen, it clearly is NOT time-invariant and there's only one value for all four cases. On top of which, it's both negative and quite a lot larger in magnitude than the computed u_s values; just from appearances it doesn't look like they'll ever get to that low a magnitude unless there's another rolloff with time.
The bottom is the ratio using the last time history value as the steady state value -- given that it is negative and the numerator positive, the ratio is negative and given the difference in magnitude, doesn't approach unity but somewhere in the neighborhood of -0.8. You'll have to figure out just how to compute a real SS value if don't want to use some variant of the empirical one. I didn't investigate whether the other issues went away with this functional; if it doesn't blow up like the other, it might.
But again, just turning off the condition warning doesn't make the numerical situation go away...
warning('off')
test_1()
function test_1()
clc; close all;
t_values = logspace(0, 2, 100);
rho1_vals = [0, 10, 50, 100]; % 4 specific curves to compare
params = struct('s1',1, 's2',100, 'c',0.000001,'beta1',1000000, 'beta2',1000000, 'b',2,'j',0.1);
figure; hold on;
hAx(1)=subplot(3,1,1); hold on % make two plots to see what are getting...
hAx(2)=subplot(3,1,2); hold on
hAx(3)=subplot(3,1,3); hold on
colors = ['r', 'g', 'b', 'm']; % Three colors for the three curves
for k = 1:length(rho1_vals)
params.rho1 = rho1_vals(k);
u_s1 = zeros(size(t_values));
u_steady = zeros(size(t_values));
for i = 1:length(t_values)
% Compute both inversions at the current beta1
u_s1(i) = talbot_inversion(@(s) U_func_s1(s, params), t_values(i));
u_steady(i) = talbot_inversion(@(s) U_func_steady(params), t_values(i));
end
% Calculate ratio
%ratio = u_s1 ./ (u_steady + eps);
ratio = u_s1/u_steady(end);
% Plot
disp([k u_s1(end) u_steady(end)])
plot(hAx(1),t_values, u_s1, 'Color', colors(k), 'LineWidth', 2, ...
'DisplayName', sprintf('\\rho=%.2f', rho1_vals(k)));
set(hAx(1), 'XScale', 'log');
xlabel(hAx(1),'Time (t)'); ylabel(hAx(1),'U_{uns}');
legend(hAx(1),'Location', 'best'); grid on;
plot(hAx(2),t_values, u_steady, 'Color', colors(k), 'LineWidth', 2, ...
'DisplayName', sprintf('\\rho=%.2f', rho1_vals(k)));
set(hAx(2), 'XScale', 'log');
xlabel(hAx(2),'Time (t)'); ylabel(hAx(2),'U_{std}');
legend(hAx(2),'Location', 'best'); grid on;
plot(hAx(3),t_values, ratio, 'Color', colors(k), 'LineWidth', 2, ...
'DisplayName', sprintf('\\rho=%.2f', rho1_vals(k)));
set(hAx(3), 'XScale', 'log');
ylim(hAx(3),[-1 0])
xlabel(hAx(3),'Time (t)'); ylabel(hAx(2),'Ratio');
legend(hAx(3),'Location', 'best'); grid on;
end
%set(gca, 'XScale', 'log');
%xlabel('Time (t)'); ylabel('Ratio (U_{uns}/U_{std})');
%legend('Location', 'best'); grid on;
% title('Effect of Micropolar Coupling on Dispersion Ratio');
end
% --- Rename your existing logic to avoid conflicts ---
function val = U_func_s1(sigma, p)
% Copy your logic from the original calc_s1.m U_function here
% Ensure it returns 'val'
% Unpack
b = p.b; c = p.c; s1 = p.s1; s2 = p.s2;
beta1 = p.beta1; beta2 = p.beta2; j = p.j; rho1 = p.rho1;
% Define S (sum) and P (product)
S = (4*c*s1 + (1 + j*s1*(c + 1))*sigma) / (c + 1);
P = (sigma*s1*(4*c + j*sigma)) / (c + 1);
% Calculate alpha1 and alpha2
alpha1 = sqrt((S - sqrt(S^2 - 4*P)) / 2);
alpha2 = sqrt((S + sqrt(S^2 - 4*P)) / 2);
%A1,B1,C1,D1,E1,F1
a11= -2 ;a12= - 2 ;
a13=- 2 *besselk(0.3e1 / 0.2e1, alpha1) ;
a14=- 2 *besselk(0.3e1 / 0.2e1, alpha2) ;
a15=- 2 *besseli(0.3e1 / 0.2e1, alpha1) ;
a16=- 2 *besseli(0.3e1 / 0.2e1, alpha2);
a21= -(beta1 + 6) / beta1 ;
a22= 2;
a23= (-alpha1 * (beta1 + 2) / beta1 * besselk(0.1e1 / 0.2e1, alpha1) - (beta1 + sigma + 6) / beta1 * besselk(0.3e1 / 0.2e1, alpha1)) ;
a24= (-alpha2 * (beta1 + 2) / beta1 * besselk(0.1e1 / 0.2e1, alpha2) - (beta1 + sigma + 6) / beta1 * besselk(0.3e1 / 0.2e1, alpha2)) ;
a25= (alpha1 * (beta1 + 2) / beta1 * besseli(0.1e1 / 0.2e1, alpha1) - (beta1 + sigma + 6) / beta1 * besseli(0.3e1 / 0.2e1, alpha1)) ;
a26= (alpha2 * (beta1 + 2) / beta1 * besseli(0.1e1 / 0.2e1, alpha2) - (beta1 + sigma + 6) / beta1 * besseli(0.3e1 / 0.2e1, alpha2)) ;
%a31=0;a32=0;
a33=((alpha1 ^ 2 * beta2 * c * s1 * s2 + alpha1 ^ 2 * beta2 * s1 * s2 - alpha1 ^ 2 * c * s1 - beta2 * s1 * s2 * sigma + alpha1 ^ 2 * c - alpha1 ^ 2 * s1 + alpha1 ^ 2 + s1 * sigma - sigma) / beta2 / s1 / c / s2 * besselk(0.3e1 / 0.2e1, alpha1) / 0.2e1 + 0.1e1 / beta2 / s2 / c * (alpha1 ^ 2 * c + alpha1 ^ 2 - sigma) * besselk(0.5e1 / 0.2e1, alpha1) * alpha1 / 0.2e1) ;
a34=((alpha2 ^ 2 * beta2 * c * s1 * s2 + alpha2 ^ 2 * beta2 * s1 * s2 - alpha2 ^ 2 * c * s1 - beta2 * s1 * s2 * sigma + alpha2 ^ 2 * c - alpha2 ^ 2 * s1 + alpha2 ^ 2 + s1 * sigma - sigma) / beta2 / s1 / c / s2 * besselk(0.3e1 / 0.2e1, alpha2) / 0.2e1 + 0.1e1 / beta2 / s2 / c * (alpha2 ^ 2 * c + alpha2 ^ 2 - sigma) * besselk(0.5e1 / 0.2e1, alpha2) * alpha2 / 0.2e1) ;
a35=((alpha1 ^ 2 * beta2 * c * s1 * s2 + alpha1 ^ 2 * beta2 * s1 * s2 - alpha1 ^ 2 * c * s1 - beta2 * s1 * s2 * sigma + alpha1 ^ 2 * c - alpha1 ^ 2 * s1 + alpha1 ^ 2 + s1 * sigma - sigma) / beta2 / s1 / c / s2 * besseli(0.3e1 / 0.2e1, alpha1) / 0.2e1 - 0.1e1 / beta2 / s2 / c * (alpha1 ^ 2 * c + alpha1 ^ 2 - sigma) * alpha1 * besseli(0.5e1 / 0.2e1, alpha1) / 0.2e1) ;
a36=((alpha2 ^ 2 * beta2 * c * s1 * s2 + alpha2 ^ 2 * beta2 * s1 * s2 - alpha2 ^ 2 * c * s1 - beta2 * s1 * s2 * sigma + alpha2 ^ 2 * c - alpha2 ^ 2 * s1 + alpha2 ^ 2 + s1 * sigma - sigma) / beta2 / s1 / c / s2 * besseli(0.3e1 / 0.2e1, alpha2) / 0.2e1 - 0.1e1 / beta2 / s2 / c * (alpha2 ^ 2 * c + alpha2 ^ 2 - sigma) * besseli(0.5e1 / 0.2e1, alpha2) * alpha2 / 0.2e1) ;
a41= -2 / b ^ 3 ;
a42=- 2 ;
a43=- 2 * b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha1) ;
a44=- 2 * b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha2) ;
a45=- 2 * b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha1) ;
a46=- 2 * b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha2);
a51= -1 / b ^ 3 ;
a52= 2 ;
a53= (-b ^ (-0.1e1 / 0.2e1) * alpha1 * besselk(0.1e1 / 0.2e1, b * alpha1) - b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha1)) ;
a54= (-b ^ (-0.1e1 / 0.2e1) * alpha2 * besselk(0.1e1 / 0.2e1, b * alpha2) - b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha2)) ;
a55= (b ^ (-0.1e1 / 0.2e1) * alpha1 * besseli(0.1e1 / 0.2e1, b * alpha1) - b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha1)) ;
a56= (b ^ (-0.1e1 / 0.2e1) * alpha2 * besseli(0.1e1 / 0.2e1, b * alpha2) - b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha2)) ;
%a61 0;a62=0;
a63=0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha1 ^ 2 * b ^ 2 * c + alpha1 ^ 2 * b ^ 2 - b ^ 2 * sigma) * besselk(0.3e1 / 0.2e1, b * alpha1) / 0.2e1 ;
a64=0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha2 ^ 2 * b ^ 2 * c + alpha2 ^ 2 * b ^ 2 - b ^ 2 * sigma) * besselk(0.3e1 / 0.2e1, b * alpha2) / 0.2e1 ;
a65=0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha1 ^ 2 * b ^ 2 * c + alpha1 ^ 2 * b ^ 2 - b ^ 2 * sigma) * besseli(0.3e1 / 0.2e1, b * alpha1)/ 0.2e1 ;
a66=0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha2 ^ 2 * b ^ 2 * c + alpha2 ^ 2 * b ^ 2 - b ^ 2 * sigma) * besseli(0.3e1 / 0.2e1, b * alpha2) / 0.2e1;
% Construct the 6x6 matrix A and column vector B
AA = [a11, a12, a13, a14, a15, a16;
a21, a22, a23, a24, a25, a26;
0, 0, a33, a34, a35, a36;
a41, a42, a43, a44, a45, a46;
a51, a52, a53, a54, a55, a56;
0, 0, a63, a64, a65, a66];
BB = [1; -1; 0; 0; 0; 0]; % Example: solve Ax=B for unit source
x = AA \ BB;
x1 = x(1);
val =( 3/(4*pi*sigma^2*(rho1 - 1 - 3*x1)));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function val0 = U_func_steady(p)
% Copy your logic from the original calc_steady.m U_function here
% Ensure it returns 'val0'
% Unpack
b = p.b; c = p.c; s1 = p.s1; s2 = p.s2;
beta1 = p.beta1; beta2 = p.beta2;
% Define S (sum) and P (product)
alpha = sqrt(4*c*s1 /(c + 1));
a11 = -2 ; a12= - 2 ;a13= - 2 ;a14= - 2 ;
a15=- 2 * besselk(0.3e1 / 0.2e1, alpha) ;
a16=- 2 * besseli(0.3e1 / 0.2e1, alpha) ;
a21 = -(beta1 + 6) / beta1 ;
a22=2 ;
a23=2 * (2 * beta1 - 5 * c - 3) / beta1 ;
a24=(beta1 + 2 * c) / beta1 ;
a25=(-alpha * (beta1 + 2) / beta1 * besselk(0.1e1 / 0.2e1, alpha) - (beta1 + 6) / beta1 * besselk(0.3e1 / 0.2e1, alpha)) ;
a26=(alpha * (beta1 + 2) / beta1 * besseli(0.1e1 / 0.2e1, alpha) - (beta1 + 6) / beta1 * besseli(0.3e1 / 0.2e1, alpha)) ;
%a31=0;a32=0;a33=0;a34=0;
a35= ((1 / beta2 * (1 + c) * alpha ^ 3 / s2 / c * besselk(0.5e1 / 0.2e1, alpha)) / 0.2e1 + ((1 + c) * (alpha ^ 2 * beta2 * s1 * s2 - alpha ^ 2 * s1 + alpha ^ 2) / beta2 / s2 / c / s1 * besselk(0.3e1 / 0.2e1, alpha)) / 0.2e1) ;
a36= (-(1 / beta2 * (1 + c) * alpha ^ 3 / s2 / c * besseli(0.5e1 / 0.2e1, alpha)) / 0.2e1 + ((1 + c) * (alpha ^ 2 * beta2 * s1 * s2 - alpha ^ 2 * s1 + alpha ^ 2) / beta2 / s2 / c / s1 * besseli(0.3e1 / 0.2e1, alpha)) / 0.2e1) ;
a41=-2 / b ^ 3 ;
a42=- 2;
a43=- 2 * b ^ 2 ;
a44=- 2 / b ;
a45=- 2 * b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha) ;
a46=- 2 * b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha);
a51= -1 / b ^ 3 ;
a52= 2 ;
a53=4 * b ^ 2 ;
a54= 1 / b ;
a55=(-b ^ (-0.1e1 / 0.2e1) * alpha * besselk(0.1e1 / 0.2e1, b * alpha) - b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha)) ;
a56=(b ^ (-0.1e1 / 0.2e1) * alpha * besseli(0.1e1 / 0.2e1, b * alpha) - b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha)) ;
%a61=0;a62=0;a63=0;a64=0;
a65= ((1 + c) / c * b ^ (-0.1e1 / 0.2e1) * alpha ^ 2 * besselk(0.3e1 / 0.2e1, b * alpha) ) / 0.2e1;
a66= ((1 + c) / c * b ^ (-0.1e1 / 0.2e1) * alpha ^ 2 * besseli(0.3e1 / 0.2e1, b * alpha)) / 0.2e1;
% Construct the 6x6 matrix A and column vector B
A0 = [a11, a12, a13, a14, a15, a16;
a21, a22, a23, a24, a25, a26;
0, 0, 0, 0, a35, a36;
a41, a42, a43, a44, a45, a46;
a51, a52, a53, a54, a55, a56;
0, 0, 0, 0, a65, a66];
B0 = [1; -1; 0; 0; 0; 0]; % Example: solve Ax=B for unit source
x0 = A0 \ B0;
x0 = x0(1);
% val0 =(3/(4*pi*sigma^2*(1 + 3*x0)));
% val0 =(3/(4*pi*(1 + 3*x0)));
val0 =1/(6*pi*x0);
end
function f = talbot_inversion(F, t)
N = 25;
f = 0;
% Talbot Contour
for k = 1:N
theta = (k-0.5)*pi/N;
s = (N/t) * (theta*cot(theta) + 1i*theta);
ds = (N/t) * (theta*(-csc(theta)^2) + cot(theta) + 1i);
% Talbot weight
f = f + exp(s*t) * F(s) * ds;
end
f = real(f / (2*pi*1i));
end
dpb
on 8 Jul 2026 at 20:20
@Torsten -- There's no real klews given as to how the end result of the "steady state" functional was generated other than comments in the third person ("Copy your logic from the original calc_steady.m U_function here") indicate that @Shreen El-Sapa got it from someone else's suggestion to do it that way, but what the differences introduced at the bottom are and how they were derived is a complete mystery.
And, of course, as the above example shows, using that functional with t generates a vector that is dependent upon the time so definitely is NOT a generic SS solution plus it is independent of the other varied parameters which one presumes it should depend upon. And, to add to the issues, as noted above, it has different sign and significantly different magnitude despite the sign than it appears any asymptotic values will be.
And in the end, since the original code divides the numerator by the time-dependent so-called but not steady-state vector, it changes the ratio with time instead of normalizing it. And, of course, that assumes one can trust the returned numeric values at all given the condition number issues. MATLAB is pretty good in dealing with that, but one shouldn't just bithely ignore the warning.
Well, I couldn't stand it...inquiring minds,and all that. Let's just see what we get if we try the empirical SS calculation with the revised functional...interestingly, it appears that they are all approaching approx 0.023 -- well, I then added 1000x and they all went up to about 0.0775...
Need to extend the evaluation time again to see what happens to the functional at larger t...
If one uses tInf = 10E5 s, that appears to be pretty good and stable range; beyond that it looks like the functional would blow up again. That's purely an empirical observation, however, dunno how one would calculate analytically.
I used 10*t(end) here for the plot; this range is 1000X the OP's original so that would be 10X my first 100X guess thinking was linear time scale...this is all empirical by observation of the behavior of the function, though. As noted, I've no idea how one would go about trying to compute independently not having any idea just what this actually is other than a bunch of numbers.
ERRATUM
The original post used the last SS value in the plot whereas the normalization used that outlined above; corrected the plot to show that set. Whether one could extend the time on out until those all did converge more closely to each other before the functional blows up with the large magnitude being passed to the Bessel functions, etc., or not I didn't test, but up to this point it does appear all do begin to converge to a similar if not identical value.
As @Torsten noted in the sidebar on the inversion code, precision issues may be a problem as well as the condition number effects on accuracy as well. It all may be just noise.
warning('off')
test_1()
function test_1()
clc; close all;
t_values = logspace(0, 4, 100);
rho1_vals = [0, 10, 50, 100]; % 4 specific curves to compare
params = struct('s1',1, 's2',100, 'c',0.000001,'beta1',1000000, 'beta2',1000000, 'b',2,'j',0.1);
figure; hold on;
hAx(1)=subplot(3,1,1); hold on % make two plots to see what are getting...
hAx(2)=subplot(3,1,2); hold on
hAx(3)=subplot(3,1,3); hold on
colors = ['r', 'g', 'b', 'm']; % Three colors for the three curves
fprintf('\n k inf realmax t(end) 10*t(end) 100*t(end) 1000*t(end)\n')
for k = 1:length(rho1_vals)
params.rho1 = rho1_vals(k);
u_s1 = zeros(size(t_values));
%u_steady = zeros(size(t_values));
for i = 1:length(t_values)
% Compute both inversions at the current beta1
u_s1(i) = talbot_inversion(@(s) U_func_s1(s, params), t_values(i));
end
u_steady(1)=talbot_inversion(@(s) U_func_s1(s,params),inf);
u_steady(2)=talbot_inversion(@(s) U_func_s1(s,params),realmax);
u_steady(3)=talbot_inversion(@(s) U_func_s1(s,params),t_values(end));
u_steady(4)=talbot_inversion(@(s) U_func_s1(s,params),10*t_values(end));
u_steady(5)=talbot_inversion(@(s) U_func_s1(s,params),100*t_values(end));
u_steady(6)=talbot_inversion(@(s) U_func_s1(s,params),1000*t_values(end));
disp([k u_steady])
%continue
% Calculate ratio
%ratio = u_s1 ./ (u_steady + eps);
%whos u_s1 u_steady ratio
SS_Val=u_steady(4); % 10E5 seconds arbitrary by inspection
ratio = u_s1/SS_Val;
% Plot
plot(hAx(1),t_values, u_s1, 'Color', colors(k), 'LineWidth', 2, ...
'DisplayName', sprintf('\\rho=%.2f', rho1_vals(k)));
set(hAx(1), 'XScale', 'log');
xlabel(hAx(1),'Time (t)'); ylabel(hAx(1),'U_{uns}');
legend(hAx(1),'Location', 'best'); grid on;
ss=repmat(SS_Val,size(t_values));
plot(hAx(2),t_values, ss, 'Color', colors(k), 'LineWidth', 2, ...
'DisplayName', sprintf('\\rho=%.2f', rho1_vals(k)));
set(hAx(2), 'XScale', 'log');
xlabel(hAx(2),'Time (t)'); ylabel(hAx(2),'U_{std}');
legend(hAx(2),'Location', 'best'); grid on;
plot(hAx(3),t_values, ratio, 'Color', colors(k), 'LineWidth', 2, ...
'DisplayName', sprintf('\\rho=%.2f', rho1_vals(k)));
set(hAx(3), 'XScale', 'log');
%ylim(hAx(3),[-1 0])
xlabel(hAx(3),'Time (t)'); ylabel(hAx(3),'Ratio');
legend(hAx(3),'Location', 'best'); grid on;
end
%set(gca, 'XScale', 'log');
%xlabel('Time (t)'); ylabel('Ratio (U_{uns}/U_{std})');
%legend('Location', 'best'); grid on;
% title('Effect of Micropolar Coupling on Dispersion Ratio');
end
% --- Rename your existing logic to avoid conflicts ---
function val = U_func_s1(sigma, p)
% Copy your logic from the original calc_s1.m U_function here
% Ensure it returns 'val'
% Unpack
b = p.b; c = p.c; s1 = p.s1; s2 = p.s2;
beta1 = p.beta1; beta2 = p.beta2; j = p.j; rho1 = p.rho1;
% Define S (sum) and P (product)
S = (4*c*s1 + (1 + j*s1*(c + 1))*sigma) / (c + 1);
P = (sigma*s1*(4*c + j*sigma)) / (c + 1);
% Calculate alpha1 and alpha2
alpha1 = sqrt((S - sqrt(S^2 - 4*P)) / 2);
alpha2 = sqrt((S + sqrt(S^2 - 4*P)) / 2);
%A1,B1,C1,D1,E1,F1
a11= -2 ;a12= - 2 ;
a13=- 2 *besselk(0.3e1 / 0.2e1, alpha1) ;
a14=- 2 *besselk(0.3e1 / 0.2e1, alpha2) ;
a15=- 2 *besseli(0.3e1 / 0.2e1, alpha1) ;
a16=- 2 *besseli(0.3e1 / 0.2e1, alpha2);
a21= -(beta1 + 6) / beta1 ;
a22= 2;
a23= (-alpha1 * (beta1 + 2) / beta1 * besselk(0.1e1 / 0.2e1, alpha1) - (beta1 + sigma + 6) / beta1 * besselk(0.3e1 / 0.2e1, alpha1)) ;
a24= (-alpha2 * (beta1 + 2) / beta1 * besselk(0.1e1 / 0.2e1, alpha2) - (beta1 + sigma + 6) / beta1 * besselk(0.3e1 / 0.2e1, alpha2)) ;
a25= (alpha1 * (beta1 + 2) / beta1 * besseli(0.1e1 / 0.2e1, alpha1) - (beta1 + sigma + 6) / beta1 * besseli(0.3e1 / 0.2e1, alpha1)) ;
a26= (alpha2 * (beta1 + 2) / beta1 * besseli(0.1e1 / 0.2e1, alpha2) - (beta1 + sigma + 6) / beta1 * besseli(0.3e1 / 0.2e1, alpha2)) ;
%a31=0;a32=0;
a33=((alpha1 ^ 2 * beta2 * c * s1 * s2 + alpha1 ^ 2 * beta2 * s1 * s2 - alpha1 ^ 2 * c * s1 - beta2 * s1 * s2 * sigma + alpha1 ^ 2 * c - alpha1 ^ 2 * s1 + alpha1 ^ 2 + s1 * sigma - sigma) / beta2 / s1 / c / s2 * besselk(0.3e1 / 0.2e1, alpha1) / 0.2e1 + 0.1e1 / beta2 / s2 / c * (alpha1 ^ 2 * c + alpha1 ^ 2 - sigma) * besselk(0.5e1 / 0.2e1, alpha1) * alpha1 / 0.2e1) ;
a34=((alpha2 ^ 2 * beta2 * c * s1 * s2 + alpha2 ^ 2 * beta2 * s1 * s2 - alpha2 ^ 2 * c * s1 - beta2 * s1 * s2 * sigma + alpha2 ^ 2 * c - alpha2 ^ 2 * s1 + alpha2 ^ 2 + s1 * sigma - sigma) / beta2 / s1 / c / s2 * besselk(0.3e1 / 0.2e1, alpha2) / 0.2e1 + 0.1e1 / beta2 / s2 / c * (alpha2 ^ 2 * c + alpha2 ^ 2 - sigma) * besselk(0.5e1 / 0.2e1, alpha2) * alpha2 / 0.2e1) ;
a35=((alpha1 ^ 2 * beta2 * c * s1 * s2 + alpha1 ^ 2 * beta2 * s1 * s2 - alpha1 ^ 2 * c * s1 - beta2 * s1 * s2 * sigma + alpha1 ^ 2 * c - alpha1 ^ 2 * s1 + alpha1 ^ 2 + s1 * sigma - sigma) / beta2 / s1 / c / s2 * besseli(0.3e1 / 0.2e1, alpha1) / 0.2e1 - 0.1e1 / beta2 / s2 / c * (alpha1 ^ 2 * c + alpha1 ^ 2 - sigma) * alpha1 * besseli(0.5e1 / 0.2e1, alpha1) / 0.2e1) ;
a36=((alpha2 ^ 2 * beta2 * c * s1 * s2 + alpha2 ^ 2 * beta2 * s1 * s2 - alpha2 ^ 2 * c * s1 - beta2 * s1 * s2 * sigma + alpha2 ^ 2 * c - alpha2 ^ 2 * s1 + alpha2 ^ 2 + s1 * sigma - sigma) / beta2 / s1 / c / s2 * besseli(0.3e1 / 0.2e1, alpha2) / 0.2e1 - 0.1e1 / beta2 / s2 / c * (alpha2 ^ 2 * c + alpha2 ^ 2 - sigma) * besseli(0.5e1 / 0.2e1, alpha2) * alpha2 / 0.2e1) ;
a41= -2 / b ^ 3 ;
a42=- 2 ;
a43=- 2 * b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha1) ;
a44=- 2 * b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha2) ;
a45=- 2 * b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha1) ;
a46=- 2 * b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha2);
a51= -1 / b ^ 3 ;
a52= 2 ;
a53= (-b ^ (-0.1e1 / 0.2e1) * alpha1 * besselk(0.1e1 / 0.2e1, b * alpha1) - b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha1)) ;
a54= (-b ^ (-0.1e1 / 0.2e1) * alpha2 * besselk(0.1e1 / 0.2e1, b * alpha2) - b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha2)) ;
a55= (b ^ (-0.1e1 / 0.2e1) * alpha1 * besseli(0.1e1 / 0.2e1, b * alpha1) - b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha1)) ;
a56= (b ^ (-0.1e1 / 0.2e1) * alpha2 * besseli(0.1e1 / 0.2e1, b * alpha2) - b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha2)) ;
%a61 0;a62=0;
a63=0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha1 ^ 2 * b ^ 2 * c + alpha1 ^ 2 * b ^ 2 - b ^ 2 * sigma) * besselk(0.3e1 / 0.2e1, b * alpha1) / 0.2e1 ;
a64=0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha2 ^ 2 * b ^ 2 * c + alpha2 ^ 2 * b ^ 2 - b ^ 2 * sigma) * besselk(0.3e1 / 0.2e1, b * alpha2) / 0.2e1 ;
a65=0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha1 ^ 2 * b ^ 2 * c + alpha1 ^ 2 * b ^ 2 - b ^ 2 * sigma) * besseli(0.3e1 / 0.2e1, b * alpha1)/ 0.2e1 ;
a66=0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha2 ^ 2 * b ^ 2 * c + alpha2 ^ 2 * b ^ 2 - b ^ 2 * sigma) * besseli(0.3e1 / 0.2e1, b * alpha2) / 0.2e1;
% Construct the 6x6 matrix A and column vector B
AA = [a11, a12, a13, a14, a15, a16;
a21, a22, a23, a24, a25, a26;
0, 0, a33, a34, a35, a36;
a41, a42, a43, a44, a45, a46;
a51, a52, a53, a54, a55, a56;
0, 0, a63, a64, a65, a66];
BB = [1; -1; 0; 0; 0; 0]; % Example: solve Ax=B for unit source
x = AA \ BB;
x1 = x(1);
val =( 3/(4*pi*sigma^2*(rho1 - 1 - 3*x1)));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function val0 = U_func_steady(p)
% Copy your logic from the original calc_steady.m U_function here
% Ensure it returns 'val0'
% Unpack
b = p.b; c = p.c; s1 = p.s1; s2 = p.s2;
beta1 = p.beta1; beta2 = p.beta2;
% Define S (sum) and P (product)
alpha = sqrt(4*c*s1 /(c + 1));
a11 = -2 ; a12= - 2 ;a13= - 2 ;a14= - 2 ;
a15=- 2 * besselk(0.3e1 / 0.2e1, alpha) ;
a16=- 2 * besseli(0.3e1 / 0.2e1, alpha) ;
a21 = -(beta1 + 6) / beta1 ;
a22=2 ;
a23=2 * (2 * beta1 - 5 * c - 3) / beta1 ;
a24=(beta1 + 2 * c) / beta1 ;
a25=(-alpha * (beta1 + 2) / beta1 * besselk(0.1e1 / 0.2e1, alpha) - (beta1 + 6) / beta1 * besselk(0.3e1 / 0.2e1, alpha)) ;
a26=(alpha * (beta1 + 2) / beta1 * besseli(0.1e1 / 0.2e1, alpha) - (beta1 + 6) / beta1 * besseli(0.3e1 / 0.2e1, alpha)) ;
%a31=0;a32=0;a33=0;a34=0;
a35= ((1 / beta2 * (1 + c) * alpha ^ 3 / s2 / c * besselk(0.5e1 / 0.2e1, alpha)) / 0.2e1 + ((1 + c) * (alpha ^ 2 * beta2 * s1 * s2 - alpha ^ 2 * s1 + alpha ^ 2) / beta2 / s2 / c / s1 * besselk(0.3e1 / 0.2e1, alpha)) / 0.2e1) ;
a36= (-(1 / beta2 * (1 + c) * alpha ^ 3 / s2 / c * besseli(0.5e1 / 0.2e1, alpha)) / 0.2e1 + ((1 + c) * (alpha ^ 2 * beta2 * s1 * s2 - alpha ^ 2 * s1 + alpha ^ 2) / beta2 / s2 / c / s1 * besseli(0.3e1 / 0.2e1, alpha)) / 0.2e1) ;
a41=-2 / b ^ 3 ;
a42=- 2;
a43=- 2 * b ^ 2 ;
a44=- 2 / b ;
a45=- 2 * b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha) ;
a46=- 2 * b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha);
a51= -1 / b ^ 3 ;
a52= 2 ;
a53=4 * b ^ 2 ;
a54= 1 / b ;
a55=(-b ^ (-0.1e1 / 0.2e1) * alpha * besselk(0.1e1 / 0.2e1, b * alpha) - b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha)) ;
a56=(b ^ (-0.1e1 / 0.2e1) * alpha * besseli(0.1e1 / 0.2e1, b * alpha) - b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha)) ;
%a61=0;a62=0;a63=0;a64=0;
a65= ((1 + c) / c * b ^ (-0.1e1 / 0.2e1) * alpha ^ 2 * besselk(0.3e1 / 0.2e1, b * alpha) ) / 0.2e1;
a66= ((1 + c) / c * b ^ (-0.1e1 / 0.2e1) * alpha ^ 2 * besseli(0.3e1 / 0.2e1, b * alpha)) / 0.2e1;
% Construct the 6x6 matrix A and column vector B
A0 = [a11, a12, a13, a14, a15, a16;
a21, a22, a23, a24, a25, a26;
0, 0, 0, 0, a35, a36;
a41, a42, a43, a44, a45, a46;
a51, a52, a53, a54, a55, a56;
0, 0, 0, 0, a65, a66];
B0 = [1; -1; 0; 0; 0; 0]; % Example: solve Ax=B for unit source
x0 = A0 \ B0;
x0 = x0(1);
% val0 =(3/(4*pi*sigma^2*(1 + 3*x0)));
% val0 =(3/(4*pi*(1 + 3*x0)));
val0 =1/(6*pi*x0);
end
function f = talbot_inversion(F, t)
N = 25;
f = 0;
% Talbot Contour
for k = 1:N
theta = (k-0.5)*pi/N;
s = (N/t) * (theta*cot(theta) + 1i*theta);
ds = (N/t) * (theta*(-csc(theta)^2) + cot(theta) + 1i);
% Talbot weight
f = f + exp(s*t) * F(s) * ds;
end
f = real(f / (2*pi*1i));
end
I'd never heard of Talbot's method for numerically computing inverse Laplace transform, so thought I'd give it a try with a simple function.
syms t real
f(t) = 1 - exp(-t);
syms s
F(s) = laplace(f(t),t,s);
tval = 0:.1:6;
for ii = 1:numel(tval)
fval(ii) = talbot_inversion(matlabFunction(F(s)),tval(ii));
end
The function fails at t = 0, which isn't suprising given the N/t terms. Is that a known limitation of Talbot inversion?
fval(1)
The rest of the points don't look too good either compared to the original signal.
figure
fplot(f(t),[0,tval(end)]);
hold on
plot(tval,fval,'o')
Is the talbot_inversion function implemented correctly?
function f = talbot_inversion(F, t)
N = 25;
f = 0;
% Talbot Contour
for k = 1:N
theta = (k-0.5)*pi/N;
s = (N/t) * (theta*cot(theta) + 1i*theta);
ds = (N/t) * (theta*(-csc(theta)^2) + cot(theta) + 1i);
% Talbot weight
f = f + exp(s*t) * F(s) * ds;
end
f = real(f / (2*pi*1i));
end
It looks to be a scaling issue, maybe?
syms t real
f(t) = 1 - exp(-t);
syms s
F(s) = laplace(f(t),t,s);
tval = 0:.1:6;
for ii = 1:numel(tval)
fval(ii) = talbot_inversion(matlabFunction(F(s)),tval(ii));
end
figure
hL=fplot(f(t),[0,tval(end)]);
hold on
plot(tval,fval,'o')
fval1=fval/(5/4*pi);
plot(tval,fval1,'x')
y=1-exp(-tval);
r=fval1./y;
plot(tval,r,'-k')
[min(r) max(r) nanmean(r) range(r)]
Well, it is almost just a linear scaling...where the almost 5/4 factor is from I don't know...just observation is what was needed.
function f = talbot_inversion(F, t)
N = 25;
f = 0;
% Talbot Contour
for k = 1:N
theta = (k-0.5)*pi/N;
s = (N/t) * (theta*cot(theta) + 1i*theta);
ds = (N/t) * (theta*(-csc(theta)^2) + cot(theta) + 1i);
% Talbot weight
f = f + exp(s*t) * F(s) * ds;
end
f = real(f / (2*pi*1i));
end
dpb
on 8 Jul 2026 at 22:26
"I'd never heard of Talbot's method ..."
I hadn't, either, but hadn't tried to find out anything but just messed with the OP's obvious difficulties in looking at what their code was doing as at least starters for him/her...
Here's what Mr. Gargel says...
"Talbot’s method is a numerical technique to invert a Laplace transform F(s) into the time domain f(t). It works by deforming the Bromwich contour into a parabola in the complex plane, which rapidly decays and avoids singularities, allowing you to use the trapezoidal rule to approximate the inverse."
Oh. That's perfectly clear -- <chuckle>
Paul
on 9 Jul 2026 at 1:06
It's quite possible that the F(s) I used in this comment does not meet some assumptions for Talbot's method as implemented in that same comment. However, I've been playing around with that Tlaobt code locally and have yet to get it to work for anything. Apparently there are different parameterizations/implementations of Talbot's method. For example, a different flavor is shown in this comment
Here is a different method @Shreen El-Sapa used before (invLaplaceStehfest). I rewrote it for high precision arithmetics (invLaplaceStehfest_hpa). Comparing results from invLaplaceStehfest and invLaplaceStehfest_hpa, I guess the Laplace inverters suffer from precision in usual arithmetic.
syms t real
f(t) = 1 - exp(-t);
syms s
F(s) = laplace(f(t),t,s);
tval = 0:.1:6;
for ii = 1:numel(tval)
fval(ii) = invLaplaceStehfest(matlabFunction(F(s)),tval(ii));
fval_hpa(ii) = invLaplaceStehfest_hpa(matlabFunction(F(s)),tval(ii));
end
figure(1)
fplot(f(t),[0,tval(end)]);
hold on
plot(tval,fval,'o')
figure(2)
fplot(f(t),[0,tval(end)]);
hold on
plot(tval,fval_hpa,'o')
% --- Stehfest Inverter ---
function f = invLaplaceStehfest(F, t, N)
if nargin < 3
N = 32;
end
ln2 = log(2); f = 0;
for n = 1:N
c_n = 0;
for k = floor((n+1)/2) : min(n, N/2)
c_n = c_n + (k^(N/2) * factorial(2*k)) / ...
(factorial(N/2 - k) * factorial(k) * factorial(k-1) * factorial(n-k) * factorial(2*k - n));
end
c_n = (-1)^(N/2 + n) * c_n;
f = f + c_n * F((n * ln2) / t);
end
f = (ln2 / t) * f;
end
% --- Stehfest Inverter ---
function f = invLaplaceStehfest_hpa(F, t, N)
if nargin < 3
N = 32;
end
digits(64);
ln2 = vpa(log(sym(2))); f = sym(0);
for nn = 1:N
n = sym(nn);
c_n = sym(0);
for kk = floor((nn+1)/2) : min(nn, N/2)
k = sym(kk);
c_n = c_n + (k^(N/2) * factorial(2*k)) / ...
(factorial(N/2 - k) * factorial(k) * factorial(k-1) * factorial(n-k) * factorial(2*k - n));
end
if mod(N/2 + nn,2) == 1
c_n = -c_n;
end
f = f + c_n * vpa(F((n * ln2) / t));
end
f = (ln2 / t) * f;
end
To prove the point about the functional eventually blowing up causing instability in the estimation of a steady-state value by direct evaluation I extended the time frame and increased the number of calculated time points to 500...then blew up the time axis to see more detail...jaggies start being visible about t=10E4 s before going completely haywire. This behavior negates the idea of using any specific multiplier or end time because don't know what the behavior is going to be. I've not looked into the root causes.
warning('off')
test_1()
hAx=gca;
figure
copyobj(hAx,gcf)
xlim([1E4 1E6])
function test_1()
t_values = logspace(0, 6, 500);
rho1_vals = [0, 10, 50, 100]; % 4 specific curves to compare
params = struct('s1',1, 's2',100, 'c',0.000001,'beta1',1000000, 'beta2',1000000, 'b',2,'j',0.1);
colors = ['r', 'g', 'b', 'm'];
for k = 1:length(rho1_vals)
params.rho1 = rho1_vals(k);
u_s1 = zeros(size(t_values));
%u_steady = zeros(size(t_values));
for i = 1:length(t_values)
% Compute both inversions at the current beta1
u_s1(i) = talbot_inversion(@(s) U_func_s1(s, params), t_values(i));
end
% Plot
semilogx(t_values, u_s1, 'Color', colors(k), 'LineWidth', 2, ...
'DisplayName', sprintf('\\rho=%.2f', rho1_vals(k)));
if k==1, hold on; end
end
xlabel('Time (t)'); ylabel('U_{uns}');
legend('Location', 'best'); grid on;
end
function val = U_func_s1(sigma, p)
% Copy your logic from the original calc_s1.m U_function here
% Ensure it returns 'val'
% Unpack
b = p.b; c = p.c; s1 = p.s1; s2 = p.s2;
beta1 = p.beta1; beta2 = p.beta2; j = p.j; rho1 = p.rho1;
% Define S (sum) and P (product)
S = (4*c*s1 + (1 + j*s1*(c + 1))*sigma) / (c + 1);
P = (sigma*s1*(4*c + j*sigma)) / (c + 1);
% Calculate alpha1 and alpha2
alpha1 = sqrt((S - sqrt(S^2 - 4*P)) / 2);
alpha2 = sqrt((S + sqrt(S^2 - 4*P)) / 2);
%A1,B1,C1,D1,E1,F1
a11= -2 ;a12= - 2 ;
a13=- 2 *besselk(0.3e1 / 0.2e1, alpha1) ;
a14=- 2 *besselk(0.3e1 / 0.2e1, alpha2) ;
a15=- 2 *besseli(0.3e1 / 0.2e1, alpha1) ;
a16=- 2 *besseli(0.3e1 / 0.2e1, alpha2);
a21= -(beta1 + 6) / beta1 ;
a22= 2;
a23= (-alpha1 * (beta1 + 2) / beta1 * besselk(0.1e1 / 0.2e1, alpha1) - (beta1 + sigma + 6) / beta1 * besselk(0.3e1 / 0.2e1, alpha1)) ;
a24= (-alpha2 * (beta1 + 2) / beta1 * besselk(0.1e1 / 0.2e1, alpha2) - (beta1 + sigma + 6) / beta1 * besselk(0.3e1 / 0.2e1, alpha2)) ;
a25= (alpha1 * (beta1 + 2) / beta1 * besseli(0.1e1 / 0.2e1, alpha1) - (beta1 + sigma + 6) / beta1 * besseli(0.3e1 / 0.2e1, alpha1)) ;
a26= (alpha2 * (beta1 + 2) / beta1 * besseli(0.1e1 / 0.2e1, alpha2) - (beta1 + sigma + 6) / beta1 * besseli(0.3e1 / 0.2e1, alpha2)) ;
%a31=0;a32=0;
a33=((alpha1 ^ 2 * beta2 * c * s1 * s2 + alpha1 ^ 2 * beta2 * s1 * s2 - alpha1 ^ 2 * c * s1 - beta2 * s1 * s2 * sigma + alpha1 ^ 2 * c - alpha1 ^ 2 * s1 + alpha1 ^ 2 + s1 * sigma - sigma) / beta2 / s1 / c / s2 * besselk(0.3e1 / 0.2e1, alpha1) / 0.2e1 + 0.1e1 / beta2 / s2 / c * (alpha1 ^ 2 * c + alpha1 ^ 2 - sigma) * besselk(0.5e1 / 0.2e1, alpha1) * alpha1 / 0.2e1) ;
a34=((alpha2 ^ 2 * beta2 * c * s1 * s2 + alpha2 ^ 2 * beta2 * s1 * s2 - alpha2 ^ 2 * c * s1 - beta2 * s1 * s2 * sigma + alpha2 ^ 2 * c - alpha2 ^ 2 * s1 + alpha2 ^ 2 + s1 * sigma - sigma) / beta2 / s1 / c / s2 * besselk(0.3e1 / 0.2e1, alpha2) / 0.2e1 + 0.1e1 / beta2 / s2 / c * (alpha2 ^ 2 * c + alpha2 ^ 2 - sigma) * besselk(0.5e1 / 0.2e1, alpha2) * alpha2 / 0.2e1) ;
a35=((alpha1 ^ 2 * beta2 * c * s1 * s2 + alpha1 ^ 2 * beta2 * s1 * s2 - alpha1 ^ 2 * c * s1 - beta2 * s1 * s2 * sigma + alpha1 ^ 2 * c - alpha1 ^ 2 * s1 + alpha1 ^ 2 + s1 * sigma - sigma) / beta2 / s1 / c / s2 * besseli(0.3e1 / 0.2e1, alpha1) / 0.2e1 - 0.1e1 / beta2 / s2 / c * (alpha1 ^ 2 * c + alpha1 ^ 2 - sigma) * alpha1 * besseli(0.5e1 / 0.2e1, alpha1) / 0.2e1) ;
a36=((alpha2 ^ 2 * beta2 * c * s1 * s2 + alpha2 ^ 2 * beta2 * s1 * s2 - alpha2 ^ 2 * c * s1 - beta2 * s1 * s2 * sigma + alpha2 ^ 2 * c - alpha2 ^ 2 * s1 + alpha2 ^ 2 + s1 * sigma - sigma) / beta2 / s1 / c / s2 * besseli(0.3e1 / 0.2e1, alpha2) / 0.2e1 - 0.1e1 / beta2 / s2 / c * (alpha2 ^ 2 * c + alpha2 ^ 2 - sigma) * besseli(0.5e1 / 0.2e1, alpha2) * alpha2 / 0.2e1) ;
a41= -2 / b ^ 3 ;
a42=- 2 ;
a43=- 2 * b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha1) ;
a44=- 2 * b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha2) ;
a45=- 2 * b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha1) ;
a46=- 2 * b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha2);
a51= -1 / b ^ 3 ;
a52= 2 ;
a53= (-b ^ (-0.1e1 / 0.2e1) * alpha1 * besselk(0.1e1 / 0.2e1, b * alpha1) - b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha1)) ;
a54= (-b ^ (-0.1e1 / 0.2e1) * alpha2 * besselk(0.1e1 / 0.2e1, b * alpha2) - b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha2)) ;
a55= (b ^ (-0.1e1 / 0.2e1) * alpha1 * besseli(0.1e1 / 0.2e1, b * alpha1) - b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha1)) ;
a56= (b ^ (-0.1e1 / 0.2e1) * alpha2 * besseli(0.1e1 / 0.2e1, b * alpha2) - b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha2)) ;
%a61 0;a62=0;
a63=0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha1 ^ 2 * b ^ 2 * c + alpha1 ^ 2 * b ^ 2 - b ^ 2 * sigma) * besselk(0.3e1 / 0.2e1, b * alpha1) / 0.2e1 ;
a64=0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha2 ^ 2 * b ^ 2 * c + alpha2 ^ 2 * b ^ 2 - b ^ 2 * sigma) * besselk(0.3e1 / 0.2e1, b * alpha2) / 0.2e1 ;
a65=0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha1 ^ 2 * b ^ 2 * c + alpha1 ^ 2 * b ^ 2 - b ^ 2 * sigma) * besseli(0.3e1 / 0.2e1, b * alpha1)/ 0.2e1 ;
a66=0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha2 ^ 2 * b ^ 2 * c + alpha2 ^ 2 * b ^ 2 - b ^ 2 * sigma) * besseli(0.3e1 / 0.2e1, b * alpha2) / 0.2e1;
% Construct the 6x6 matrix A and column vector B
AA = [a11, a12, a13, a14, a15, a16;
a21, a22, a23, a24, a25, a26;
0, 0, a33, a34, a35, a36;
a41, a42, a43, a44, a45, a46;
a51, a52, a53, a54, a55, a56;
0, 0, a63, a64, a65, a66];
BB = [1; -1; 0; 0; 0; 0]; % Example: solve Ax=B for unit source
x = AA \ BB;
x1 = x(1);
val =( 3/(4*pi*sigma^2*(rho1 - 1 - 3*x1)));
end
function f = talbot_inversion(F, t)
N = 25;
f = 0;
% Talbot Contour
for k = 1:N
theta = (k-0.5)*pi/N;
s = (N/t) * (theta*cot(theta) + 1i*theta);
ds = (N/t) * (theta*(-csc(theta)^2) + cot(theta) + 1i);
% Talbot weight
f = f + exp(s*t) * F(s) * ds;
end
f = real(f / (2*pi*1i));
end
I haven't been able to find a really good reference that explains Talbot's method that easily translates to code. The one that I did find is this ACM paper. I was able to implement the equations in that paper, but changed the second equation in (2.5) that didn't seem correct to me. Also, I tried implementing the integral in the first equation in (2.5) (using with what I think is the correct expression for s_v_prime) with integral but was not successful. I'd be very interested if someone could get that to work.
From my very, very limited understanding of the problem, the shape parameters (nu, sigma, lambda) have to be chosen to meet certain criteria based on the singularities of F(s), i.e., a single set of parameters won't necessarily work for all problems. I think that one such criterion is that all poles of F(s) be to the left of the Talbot contour s(theta;nu,sigma,lambda). Also, the approximate inversion becomes more problematic as t increases, which should be of interest in this problem that is concerned with steady state behavior.
For example, consider the difference of two decaying exponentials
F = @(s) 1./(s + 2) - 1./(s + 1);
Then we have
tval = .1:.1:6;
f0 = talbot0(F,tval,25,1,.5,25./tval);
figure
plot(tval,exp(-2*tval)-exp(-1*tval),tval,real(f0),'o'),grid
The inversion as t gets larger can be improved my reducing sigma, yet keeping it to the right of the rightmost pole of F(s)
f0 = talbot0(F,tval,25,1,-0.5,25./tval);
figure
plot(tval,exp(-2*tval)-exp(-1*tval),tval,real(f0),'o'),grid
Similarly, consider the case where f(t) is non-zero in steady state
F= @(s) 1./s - 1./(s + 1);
f0 = talbot0(F,tval,25,1,1,25./tval);
figure
plot(tval,1-exp(-1*tval),tval,real(f0),'o'),grid
I don't know why the inversion diverges for these parameters. It can be improved by reducing sigma, but keeping it to the right of the rightmost pole
f0 = talbot0(F,tval,25,1,0.1,25./tval);
figure
plot(tval,1 - exp(-1*tval),tval,real(f0),'o'),grid
It would not suprise me if the inversion diverges for larger values of t for this set of parameters.
Given these complications even for such simple problems, it seems it could be rather difficult to use this method for an F(s) that's only known numerically without having much analytical insight into F(s) or have some way to extract properties of F(s) numerically to select good parameters, as I think is done in the linked paper.
function f = talbot0(F,t,N,nu,sigma,lambda)
k = (1:N-1).';
theta = k*pi/N;
sv = theta.*cot(theta) + 1i*nu*theta;
%svp = 1i*(nu + (theta - cos(theta).*sin(theta))./sin(theta).^2);
svp = 1i*nu + cot(theta) - theta.*csc(theta).^2;
Tn = sum(exp(lambda.*t.*sv).*F(sigma + lambda.*sv).*svp./1i,1);
Tn = Tn + nu/2*exp(lambda.*t).*F(sigma+lambda);
f = lambda.*exp(sigma*t)/N.*Tn;
end
Torsten
on 10 Jul 2026 at 14:12
This code is from the File Exchange:
F = @(s) 1./(s + 2) - 1./(s + 1);
tval = .1:.1:6;
f0 = talbot_inversion(F,tval);
figure
plot(tval,exp(-2*tval)-exp(-1*tval),tval,real(f0),'o'),grid
F= @(s) 1./s - 1./(s + 1);
f0 = talbot_inversion(F,tval);
figure
plot(tval,1-exp(-1*tval),tval,real(f0),'o'),grid
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.
%
% Abate, 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
As the FEX function stated in the comments, the parameter M can be important.
Consider
wn = 1;
F = @(s) wn^2./s./(s.^2 + 2.*0.1*wn.*s + wn^2);
tval = .1:.5:30;
f0 = talbot0(F,tval,25,1,0,25./tval);
filt = talbot_inversion(F,tval,25); % Use M = 25, but works fine with default M = 64;
figure
step(tf(wn^2,[1,2*.1*wn,wn^2]),tval(end));
hold on
plot(tval,real(f0),'rx');
plot(tval,filt,'go')
grid
function f = talbot0(F,t,N,nu,sigma,lambda)
k = (1:N-1).';
theta = k*pi/N;
sv = theta.*cot(theta) + 1i*nu*theta;
%svp = 1i*(nu + (theta - cos(theta).*sin(theta))./sin(theta).^2);
svp = 1i*nu + cot(theta) - theta.*csc(theta).^2;
Tn = sum(exp(lambda.*t.*sv).*F(sigma + lambda.*sv).*svp./1i,1);
Tn = Tn + nu/2*exp(lambda.*t).*F(sigma+lambda);
f = lambda.*exp(sigma*t)/N.*Tn;
end
function ilt = talbot_inversion(f_s, t, M)
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
David Goodmanson
on 11 Jul 2026 at 10:56
Hi Paul,
you were looking to see the integral function work for the Talbot method, so here is an example. I used your exponential difference function but reversed the ovrall sign so that the semilog plot looks better. I just took some guesses for the parameters. It appears that the Talbot method was done for the purpose of quadratures but nowdays it seems like there are easier ways to make a path and just sum over a lot of points, like maybe millions.
t = .01:.1:16;
sigma = 1;
lambda = 1;
nu = 1/2;
A = zeros(size(t));
for k = 1:length(t)
tk = t(k);
I = integral(@(th) fun(th,tk,sigma,lambda,nu),-pi+.001,pi-.001);
A(k) = lambda*exp(sigma*tk)*I/(2*pi*i);
end
figure(1)
subplot(1,2,1)
plot(t,exp(-t)-exp(-2*t),t,real(A),'o')
subplot(1,2,2)
semilogy(t,exp(-t)-exp(-2*t),t,real(A),'o')
function g = fun(th,t,sigma,lambda,nu)
sbar = th./tan(th)+i*nu*th;
sbar(th==0) = 1; % this safety play doesn't seem to make any difference
ss = sigma + lambda*sbar;
dsbar = (i*nu + 1./tan(th) - th./sin(th).^2);
dsbar(th==0) = (i*nu); % neither does this one
g = exp(lambda*t*sbar).*dsbar.*(1./(ss+1) - 1./(ss+2));
end
Hi David,
Thanks for posting that solution. It inspired me to go back and check my own code and, lo and behold, it works fine for the example at hand, albeit with lots of warnings. I'm not sure why I thought differently when I last tried it.
F = @(s) 1./(s + 2) - 1./(s + 1);
tval = .1:.1:6;
warning('off')
f2 = talbot2(F,tval,25,1,0,25./tval);
figure
plot(tval,exp(-2*tval) - exp(-tval),tval,f2,'o'),grid
function f = talbot2(F,t,N,nu,sigma,lambda)
f = nan(size(t));
sv = @(theta) theta.*(cot(theta) + 1i*nu);
svp = @(theta) 1i*nu + cot(theta) - theta.*csc(theta).^2;
for ii = 1:numel(t)
I = @(theta) real(lambda(ii)/2/pi/1i*exp(sigma*t(ii))*exp(lambda(ii)*sv(theta)*t(ii)).*F(sigma + lambda(ii)*sv(theta)).*svp(theta));
f(ii) = integral(I,-pi+eps,0-eps) + integral(I,0+eps,pi-eps);
end
end
If we know the pole locations (and mabye F(s) satisfies some other criteria?), then we can evaluate the ILT integral directly by specifying an appropriate contour using the Waypoints option of integral
F = @(s) 1./(s+2) - 1./(s+1);
tval = 0.1:0.1:6;
for ii = 1:numel(tval)
f2(ii) = real(integral(@(s) exp(s*tval(ii)).*F(s),0,0,'Waypoints',[1i,-5+1i,-5-1i,-1i])/2/pi/1i);
end
figure
plot(tval,exp(-2*tval)-exp(-tval),tval,f2,'-o'),grid
Of course we don't know the pole locations for the original problem in this question, but I suspect that "knowing the poles" is conceptually related to determining the appropriate shape parameters to specify the Talbot contour in Talbot's method.
p.s. No idea why the solid line in the plot isn't blue.
David Goodmanson
on 11 Jul 2026 at 21:22
From the paper that you cited:
" A prescription for finding a suitable contour (one which retains a proper balance between the two effects described above) has been provided by Talbot. This requires that the locations of the singularities of F(s) be known, "
along with some other conditions. So Talbot is no different.
If that's truly the case, what makes the Talbot contour useful? Is it just that Talbot's contour produces a delta_s that increases as s ->> -inf, so that summation out to -inf is practical? Having a function like integral that chooses delta_s dynamically would appear to make Talbot unnecessary, although F(s) with branch cuts means you have to pay attention to the falloff of F(s) and pick an appropriate path.
p.s. Everyone here including Talbot is assuming that F(s) is analytic.
Paul
on 12 Jul 2026 at 0:09
No idea if Talbot's approach is useful. I just learned about it from this thread.
Yes, F(s) has to be analytic, and I suspect that F(s) also has to decay to zero (uniformly?) as s -> infinity.
Question: Given a software function that takes in complex s and returns a numerical value of F(s), is there a way to determine if F(s) is analytic?
Categories
Find more on Particle & Nuclear Physics in Help Center and File Exchange
on 8 Jul 2026 at 11:49
on 12 Jul 2026 at 0:09
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
