inversion Laplace to val function
Mostra commenti meno recenti
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?
5 Commenti
Shreen El-Sapa
5 minuti fa
The axis is negative because the calculation (whatever it is) returns negative values; we don't know anything about that == and there's absolutely no context provided; how are we supposed to know anything whatsoever about the problem this code is supposed to be solving?
What parameters are those? There's a list
% 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;
identified although we have no idea what they might represent. It would seem putting the start, stop times desired and the params struct in the function argument list instead of hardcoding them inside the function would be more user amenable to making changes.
But, what those are and might could be, we've no klew and no way to know...although it looks somewhat like an exponential decay process, maybe???
From the questions, one infers this is not code you wrote; the first step would be to go to the source from whence it came.
Here is a way on how you can vary your parameters.
Of course, it's impossible for us to interprete the results of your computations and thus to answer your question about the negative y-values.
clc;
close all;
t_values = linspace(0.0, 10, 100);
u_t = zeros(size(t_values));
n_params = 4;
Params.c = [1e-10;0.5;1;10];
Params.b = [10;10;10;10];
Params.s1 = [3;3;3;3];
Params.s2 = [3;3;3;3];
Params.beta1 = [1066;1066;1066;1066];
Params.beta2 = [1e6;1e6;1e6;1e6];
Params.j = [0.1;0.1;0.1;0.1];
Params.rho1 = [0;0;0;0];
hold on
for k = 1:n_params
% Define parameters that are needed in U_function
params.c = Params.c(k);
params.b = Params.b(k);
params.s1 = Params.s1(k);
params.s2 = Params.s2(k);
params.beta1 = Params.beta1(k);
params.beta2 = Params.beta2(k);
params.j = Params.j(k);
params.rho1 = Params.rho1(k);
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);
end
hold off
grid on; xlabel('Time (t)'); ylabel('u(t)');
title('Inverse Laplace Transform Results (Talbot)');
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
Shreen El-Sapa
circa 6 ore fa
One way to make the result look positive is to plot abs(u_t) instead of real(u_t)...of course, it then looks like a growth phenomenon instead of decay...
"How can I increase values of time on x-axis?"
Just change the limits in the linspace call to set t_values vector...
Tend=20; % define an end time variable instead of fixed number
t_values = linspace(0.0, Tend, 100); % use it instead
u_t = zeros(size(t_values));
n_params = 4;
Params.c = [1e-10;0.5;1;10];
Params.b = [10;10;10;10];
Params.s1 = [3;3;3;3];
Params.s2 = [3;3;3;3];
Params.beta1 = [1066;1066;1066;1066];
Params.beta2 = [1e6;1e6;1e6;1e6];
Params.j = [0.1;0.1;0.1;0.1];
Params.rho1 = [0;0;0;0];
hold on
for k = 1:n_params
% Define parameters that are needed in U_function
params.c = Params.c(k);
params.b = Params.b(k);
params.s1 = Params.s1(k);
params.s2 = Params.s2(k);
params.beta1 = Params.beta1(k);
params.beta2 = Params.beta2(k);
params.j = Params.j(k);
params.rho1 = Params.rho1(k);
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);
plot(t_values, abs(u_t), 'LineWidth', 2);
end
hold off
grid on; xlabel('Time (t)'); ylabel('u(t)');
title('Inverse Laplace Transform Results (Talbot)');
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
Risposte (0)
circa 11 ore fa
circa 5 ore fa
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
