test0
function test0()
clc;
close all;
t_values = linspace(0.0, 10, 100);
u_t = zeros(size(t_values));
% Define parameters that are needed in U_function
params.c = 10^(-10);
params.b = 10;
params.s1 = 3;
params.s2 = 3;
params.beta1 = 1066;
params.beta2 = 10^6;
params.j = 0.1;
params.rho1 = 0;
for i = 1:length(t_values)
u_t(i) = talbot_inversion(@(s) U_function(s, params), t_values(i));
end
plot(t_values, real(u_t), ‘LineWidth’, 2);
grid on; xlabel(‘Time (t)’); ylabel(‘u(t)’);
title(‘Inverse Laplace Transform Results (Talbot)’);
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.1e1 / 0.2e1, alpha1) – (4 * beta1 * c + beta1 * sigma + 2 * beta1 + 1) * besselk(0.3e1 / 0.2e1, alpha1)) ;
a24= (-alpha2 * (2 * beta1 + 1) * besselk(0.1e1 / 0.2e1, alpha2) – (4 * beta1 * c + beta1 * sigma + 2 * beta1 + 1) * besselk(0.3e1 / 0.2e1, alpha2)) ;
a25= (alpha1 * (2 * beta1 + 1) * besseli(0.1e1 / 0.2e1, alpha1) – (4 * beta1 * c + beta1 * sigma + 2 * beta1 + 1) * besseli(0.3e1 / 0.2e1, alpha1));
a26= (alpha2 * (2 * beta1 + 1) * besseli(0.1e1 / 0.2e1, alpha2) – (4 * beta1 * c + beta1 * sigma + 2 * beta1 + 1) * besseli(0.3e1 / 0.2e1, 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(0.3e1 / 0.2e1, alpha1) / 0.2e1 + beta2 / s1 / c * (alpha1 ^ 2 * c + alpha1 ^ 2 – sigma) * besselk(0.5e1 / 0.2e1, alpha1) * alpha1 / 0.2e1) ;
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(0.3e1 / 0.2e1, alpha2) / 0.2e1 + beta2 / s1 / c * (alpha2 ^ 2 * c + alpha2 ^ 2 – sigma) * besselk(0.5e1 / 0.2e1, alpha2) * alpha2 / 0.2e1) ;
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(0.3e1 / 0.2e1, alpha1) / 0.2e1 – beta2 / s1 / c * (alpha1 ^ 2 * c + alpha1 ^ 2 – sigma) * alpha1 * besseli(0.5e1 / 0.2e1, alpha1) / 0.2e1) ;
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(0.3e1 / 0.2e1, alpha2) / 0.2e1 – beta2 / s1 / 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^(-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)
%num_1 = (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);
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
So, what’s the question?test0
function test0()
clc;
close all;
t_values = linspace(0.0, 10, 100);
u_t = zeros(size(t_values));
% Define parameters that are needed in U_function
params.c = 10^(-10);
params.b = 10;
params.s1 = 3;
params.s2 = 3;
params.beta1 = 1066;
params.beta2 = 10^6;
params.j = 0.1;
params.rho1 = 0;
for i = 1:length(t_values)
u_t(i) = talbot_inversion(@(s) U_function(s, params), t_values(i));
end
plot(t_values, real(u_t), ‘LineWidth’, 2);
grid on; xlabel(‘Time (t)’); ylabel(‘u(t)’);
title(‘Inverse Laplace Transform Results (Talbot)’);
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.1e1 / 0.2e1, alpha1) – (4 * beta1 * c + beta1 * sigma + 2 * beta1 + 1) * besselk(0.3e1 / 0.2e1, alpha1)) ;
a24= (-alpha2 * (2 * beta1 + 1) * besselk(0.1e1 / 0.2e1, alpha2) – (4 * beta1 * c + beta1 * sigma + 2 * beta1 + 1) * besselk(0.3e1 / 0.2e1, alpha2)) ;
a25= (alpha1 * (2 * beta1 + 1) * besseli(0.1e1 / 0.2e1, alpha1) – (4 * beta1 * c + beta1 * sigma + 2 * beta1 + 1) * besseli(0.3e1 / 0.2e1, alpha1));
a26= (alpha2 * (2 * beta1 + 1) * besseli(0.1e1 / 0.2e1, alpha2) – (4 * beta1 * c + beta1 * sigma + 2 * beta1 + 1) * besseli(0.3e1 / 0.2e1, 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(0.3e1 / 0.2e1, alpha1) / 0.2e1 + beta2 / s1 / c * (alpha1 ^ 2 * c + alpha1 ^ 2 – sigma) * besselk(0.5e1 / 0.2e1, alpha1) * alpha1 / 0.2e1) ;
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(0.3e1 / 0.2e1, alpha2) / 0.2e1 + beta2 / s1 / c * (alpha2 ^ 2 * c + alpha2 ^ 2 – sigma) * besselk(0.5e1 / 0.2e1, alpha2) * alpha2 / 0.2e1) ;
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(0.3e1 / 0.2e1, alpha1) / 0.2e1 – beta2 / s1 / c * (alpha1 ^ 2 * c + alpha1 ^ 2 – sigma) * alpha1 * besseli(0.5e1 / 0.2e1, alpha1) / 0.2e1) ;
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(0.3e1 / 0.2e1, alpha2) / 0.2e1 – beta2 / s1 / 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^(-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)
%num_1 = (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);
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
So, what’s the question? test0
function test0()
clc;
close all;
t_values = linspace(0.0, 10, 100);
u_t = zeros(size(t_values));
% Define parameters that are needed in U_function
params.c = 10^(-10);
params.b = 10;
params.s1 = 3;
params.s2 = 3;
params.beta1 = 1066;
params.beta2 = 10^6;
params.j = 0.1;
params.rho1 = 0;
for i = 1:length(t_values)
u_t(i) = talbot_inversion(@(s) U_function(s, params), t_values(i));
end
plot(t_values, real(u_t), ‘LineWidth’, 2);
grid on; xlabel(‘Time (t)’); ylabel(‘u(t)’);
title(‘Inverse Laplace Transform Results (Talbot)’);
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.1e1 / 0.2e1, alpha1) – (4 * beta1 * c + beta1 * sigma + 2 * beta1 + 1) * besselk(0.3e1 / 0.2e1, alpha1)) ;
a24= (-alpha2 * (2 * beta1 + 1) * besselk(0.1e1 / 0.2e1, alpha2) – (4 * beta1 * c + beta1 * sigma + 2 * beta1 + 1) * besselk(0.3e1 / 0.2e1, alpha2)) ;
a25= (alpha1 * (2 * beta1 + 1) * besseli(0.1e1 / 0.2e1, alpha1) – (4 * beta1 * c + beta1 * sigma + 2 * beta1 + 1) * besseli(0.3e1 / 0.2e1, alpha1));
a26= (alpha2 * (2 * beta1 + 1) * besseli(0.1e1 / 0.2e1, alpha2) – (4 * beta1 * c + beta1 * sigma + 2 * beta1 + 1) * besseli(0.3e1 / 0.2e1, 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(0.3e1 / 0.2e1, alpha1) / 0.2e1 + beta2 / s1 / c * (alpha1 ^ 2 * c + alpha1 ^ 2 – sigma) * besselk(0.5e1 / 0.2e1, alpha1) * alpha1 / 0.2e1) ;
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(0.3e1 / 0.2e1, alpha2) / 0.2e1 + beta2 / s1 / c * (alpha2 ^ 2 * c + alpha2 ^ 2 – sigma) * besselk(0.5e1 / 0.2e1, alpha2) * alpha2 / 0.2e1) ;
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(0.3e1 / 0.2e1, alpha1) / 0.2e1 – beta2 / s1 / c * (alpha1 ^ 2 * c + alpha1 ^ 2 – sigma) * alpha1 * besseli(0.5e1 / 0.2e1, alpha1) / 0.2e1) ;
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(0.3e1 / 0.2e1, alpha2) / 0.2e1 – beta2 / s1 / 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^(-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)
%num_1 = (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);
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
So, what’s the question? laplace, inversion, function MATLAB Answers — New Questions
