plot after taking inversion Laplace to val function and then divide it on val0 function for steady case, it should give 1 as t tends to infty
function main_calculation()
clc; close all;
t_values = logspace(0, 2, 100);
beta1_vals = [0.001, 1, 3, 50]; % 4 specific curves to compare
params = struct(‘c’,0.0000001, ‘s1’,0.00001, ‘s2’,0.1, ‘rho1’,0, ‘beta2’,10000, ‘b’,2,’j’,0.1);
figure; hold on;
colors = [‘r’, ‘g’, ‘b’, ‘m’]; % Three colors for the three curves
for k = 1:length(beta1_vals)
params.beta1 = beta1_vals(k);
u_s1 = zeros(size(t_values));
u_steady = zeros(size(t_values));
for i = 1:length(t_values)
% Compute both inversions at the current beta1
u_s1(i) = talbot_inversion(@(s) U_func_s1(s, params), t_values(i));
u_steady(i) = talbot_inversion(@(s) U_func_steady(s, params), t_values(i));
end
% Calculate ratio
ratio = u_s1 ./ (u_steady + eps);
% Plot
plot(t_values, ratio, ‘Color’, colors(k), ‘LineWidth’, 2, …
‘DisplayName’, sprintf(‘\beta_1=%.2f’, beta1_vals(k)));
end
set(gca, ‘XScale’, ‘log’);
xlabel(‘Time (t)’); ylabel(‘Ratio (U_{osc}/U_{std})’);
legend(‘Location’, ‘best’); grid on;
% title(‘Effect of Micropolar Coupling on Dispersion Ratio’);
end
% — Rename your existing logic to avoid conflicts —
function val = U_func_s1(sigma, p)
% Copy your logic from the original calc_s1.m U_function here
% Ensure it returns ‘val’
% Unpack
b = p.b; c = p.c; s1 = p.s1; s2 = p.s2;
beta1 = p.beta1; beta2 = p.beta2; j = p.j; rho1 = p.rho1;
% Define S (sum) and P (product)
S = (4*c*s1 + (1 + j*s1*(c + 1))*sigma) / (c + 1);
P = (sigma*s1*(4*c + j*sigma)) / (c + 1);
% Calculate alpha1 and alpha2
alpha1 = sqrt((S – sqrt(S^2 – 4*P)) / 2);
alpha2 = sqrt((S + sqrt(S^2 – 4*P)) / 2);
%A1,B1,C1,D1,E1,F1
a11= 2; a12= 2 ;
a13= 2*besselk(0.3e1 / 0.2e1, alpha1);a14= 2 * besselk(0.3e1 / 0.2e1, alpha2) ;
a15= 2*besseli(0.3e1 / 0.2e1, alpha1);a16= 2 * besseli(0.3e1 / 0.2e1, alpha2);
a21= -(beta1 + 4 * c + 2) / beta1 ;
a22= 2 * (beta1 – 2 * c + 2) / beta1 ;
a23=(-alpha1 * (beta1 + 2) / beta1 * besselk(0.1e1 / 0.2e1, alpha1) – (beta1 + 4 * c + sigma + 2) / beta1 * besselk(0.3e1 / 0.2e1, alpha1));
a24=(-alpha2 * (beta1 + 2) / beta1 * besselk(0.1e1 / 0.2e1, alpha2) – (beta1 + 4 * c + sigma + 2) / beta1 * besselk(0.3e1 / 0.2e1, alpha2));
a25=(alpha1 * (beta1 + 2) / beta1 * besseli(0.1e1 / 0.2e1, alpha1) – (beta1 + 4 * c + sigma + 2) / beta1 * besseli(0.3e1 / 0.2e1, alpha1)) ;
a26=(alpha2 * (beta1 + 2) / beta1 * besseli(0.1e1 / 0.2e1, alpha2) – (beta1 + 4 * c + sigma + 2) / beta1 * besseli(0.3e1 / 0.2e1, alpha2));
%a31=0;a32=0;
a33= ((alpha1 ^ 2 * beta2 * c * s1 * s2 + alpha1 ^ 2 * beta2 * s1 * s2 – alpha1 ^ 2 * c * s1 – beta2 * s1 * s2 * sigma + alpha1 ^ 2 * c – alpha1 ^ 2 * s1 + alpha1 ^ 2 + s1 * sigma – sigma) / beta2 / s1 / c / s2 * besselk(0.3e1 / 0.2e1, alpha1) / 0.2e1 + 0.1e1 / beta2 / s2 / c * (alpha1 ^ 2 * c + alpha1 ^ 2 – sigma) * besselk(0.5e1 / 0.2e1, alpha1) * alpha1 / 0.2e1) ;
a34= ((alpha2 ^ 2 * beta2 * c * s1 * s2 + alpha2 ^ 2 * beta2 * s1 * s2 – alpha2 ^ 2 * c * s1 – beta2 * s1 * s2 * sigma + alpha2 ^ 2 * c – alpha2 ^ 2 * s1 + alpha2 ^ 2 + s1 * sigma – sigma) / beta2 / s1 / c / s2 * besselk(0.3e1 / 0.2e1, alpha2) / 0.2e1 + 0.1e1 / beta2 / s2 / c * (alpha2 ^ 2 * c + alpha2 ^ 2 – sigma) * besselk(0.5e1 / 0.2e1, alpha2) * alpha2 / 0.2e1) ;
a35= ((alpha1 ^ 2 * beta2 * c * s1 * s2 + alpha1 ^ 2 * beta2 * s1 * s2 – alpha1 ^ 2 * c * s1 – beta2 * s1 * s2 * sigma + alpha1 ^ 2 * c – alpha1 ^ 2 * s1 + alpha1 ^ 2 + s1 * sigma – sigma) / beta2 / s1 / c / s2 * besseli(0.3e1 / 0.2e1, alpha1) / 0.2e1 – 0.1e1 / beta2 / s2 / c * (alpha1 ^ 2 * c + alpha1 ^ 2 – sigma) * alpha1 * besseli(0.5e1 / 0.2e1, alpha1) / 0.2e1) ;
a36= ((alpha2 ^ 2 * beta2 * c * s1 * s2 + alpha2 ^ 2 * beta2 * s1 * s2 – alpha2 ^ 2 * c * s1 – beta2 * s1 * s2 * sigma + alpha2 ^ 2 * c – alpha2 ^ 2 * s1 + alpha2 ^ 2 + s1 * sigma – sigma) / beta2 / s1 / c / s2 * besseli(0.3e1 / 0.2e1, alpha2) / 0.2e1 – 0.1e1 / beta2 / s2 / c * (alpha2 ^ 2 * c + alpha2 ^ 2 – sigma) * besseli(0.5e1 / 0.2e1, alpha2) * alpha2 / 0.2e1) ;
a41=2 / b ^ 3; a42= 2;
a43= 2 * b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha1) ;
a44= 2 * b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha2);
a45= 2 * b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha1) ;
a46= 2 * b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha2);
a51= -1 / b ^ 3 ;a52= 2 ;
a53= (-b ^ (-0.1e1 / 0.2e1) * alpha1 * besselk(0.1e1 / 0.2e1, b * alpha1) – b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha1)) ;
a54= (-b ^ (-0.1e1 / 0.2e1) * alpha2 * besselk(0.1e1 / 0.2e1, b * alpha2) – b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha2)) ;
a55= (b ^ (-0.1e1 / 0.2e1) * alpha1 * besseli(0.1e1 / 0.2e1, b * alpha1) – b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha1)) ;
a56= (b ^ (-0.1e1 / 0.2e1) * alpha2 * besseli(0.1e1 / 0.2e1, b * alpha2) – b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha2)) ;
%a61=0; a62=0;
a63= 0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha1 ^ 2 * b ^ 2 * c + alpha1 ^ 2 * b ^ 2 – b ^ 2 * sigma) * besselk(0.3e1 / 0.2e1, b * alpha1)/ 0.2e1;
a64= 0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha2 ^ 2 * b ^ 2 * c + alpha2 ^ 2 * b ^ 2 – b ^ 2 * sigma) * besselk(0.3e1 / 0.2e1, b * alpha2) / 0.2e1;
a65= 0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha1 ^ 2 * b ^ 2 * c + alpha1 ^ 2 * b ^ 2 – b ^ 2 * sigma) * besseli(0.3e1 / 0.2e1, b * alpha1) / 0.2e1 ;
a66= 0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha2 ^ 2 * b ^ 2 * c + alpha2 ^ 2 * b ^ 2 – b ^ 2 * sigma) * besseli(0.3e1 / 0.2e1, b * alpha2) / 0.2e1;
% Construct the 6×6 matrix A and column vector B
A = [a11, a12, a13, a14, a15, a16;
a21, a22, a23, a24, a25, a26;
0, 0, a33, a34, a35, a36;
a41, a42, a43, a44, a45, a46;
a51, a52, a53, a54, a55, a56;
0, 0, a63, a64, a65, a66];
B = [1; -1; 0; 0; 0; 0]; % Example: solve Ax=B for unit source
x = A B;
x1 = x(1);
val =( 3/(4*pi*sigma^2*(rho1 – 1 – 3*x1)));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function val0 = U_func_steady(sigma, p)
% Copy your logic from the original calc_steady.m U_function here
% Ensure it returns ‘val0’
% Unpack
b = p.b; c = p.c; s1 = p.s1; s2 = p.s2;
beta1 = p.beta1; beta2 = p.beta2;
% Define S (sum) and P (product)
alpha = sqrt(4*c*s1 /(c + 1));
a11= 2; a12= 2 ; a13= 2 ; a14=2; a15=2 * besselk(0.3e1 / 0.2e1, alpha); a16=2 * besseli(0.3e1 / 0.2e1, alpha);
a21 = -(beta1 + 4 * c + 2) / beta1 ;
a22=2 * (beta1 – 2 * c + 2) / beta1 ;
a23=2 * (2 * beta1 – 7 * c – 1) / beta1;
a24=(beta1 – 2 * c + 4) / beta1 ;
a25=(-alpha * (beta1 + 2) / beta1 * besselk(0.1e1 / 0.2e1, alpha) – (beta1 + 4 * c + sigma + 2) / beta1 * besselk(0.3e1 / 0.2e1, alpha));
a26=(alpha * (beta1 + 2) / beta1 * besseli(0.1e1 / 0.2e1, alpha) – (beta1 + 4 * c + sigma + 2) / beta1 * besseli(0.3e1 / 0.2e1, alpha));
%a31=0;a32=0;a33=0;a34=0;
a35= ((alpha ^ 2 * beta2 * c * s1 * s2 + alpha ^ 2 * beta2 * s1 * s2 – alpha ^ 2 * c * s1 – beta2 * s1 * s2 * sigma + alpha ^ 2 * c – alpha ^ 2 * s1 + alpha ^ 2 + s1 * sigma – sigma) / beta2 / s1 / c / s2 * besselk(0.3e1 / 0.2e1, alpha) / 0.2e1 + 0.1e1 / beta2 / s2 / c * (alpha ^ 2 * c + alpha ^ 2 – sigma) * besselk(0.5e1 / 0.2e1, alpha) * alpha / 0.2e1) ;
a36= ((alpha ^ 2 * beta2 * c * s1 * s2 + alpha ^ 2 * beta2 * s1 * s2 – alpha ^ 2 * c * s1 – beta2 * s1 * s2 * sigma + alpha ^ 2 * c – alpha ^ 2 * s1 + alpha ^ 2 + s1 * sigma – sigma) / beta2 / s1 / c / s2 * besseli(0.3e1 / 0.2e1, alpha) / 0.2e1 – 0.1e1 / beta2 / s2 / c * (alpha ^ 2 * c + alpha ^ 2 – sigma) * besseli(0.5e1 / 0.2e1, alpha) * alpha / 0.2e1) ;
a41=2 / b ^ 3; a42=2 ; a43=2 * b ^ 2 ; a44=2 / b ;
a45=2 * b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha);
a46=2 * b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha);
a51= -1 / b ^ 3; a52= 2 ; a53=4 * b ^ 2 ;a54=1 / b ;
a55=(-b ^ (-0.1e1 / 0.2e1) * alpha * besselk(0.1e1 / 0.2e1, b * alpha) – b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha)) ;
a56=(b ^ (-0.1e1 / 0.2e1) * alpha * besseli(0.1e1 / 0.2e1, b * alpha) – b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha)) ;
%a61=0;a62=0;a63=0;a64=0;
a65= 0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha ^ 2 * b ^ 2 * c + alpha ^ 2 * b ^ 2 – b ^ 2 * sigma) * besselk(0.3e1 / 0.2e1, b * alpha) / 0.2e1;
a66= 0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha ^ 2 * b ^ 2 * c + alpha ^ 2 * b ^ 2 – b ^ 2 * sigma) * besseli(0.3e1 / 0.2e1, b * alpha) / 0.2e1;
% Construct the 6×6 matrix A and column vector B
A0 = [a11, a12, a13, a14, a15, a16;
a21, a22, a23, a24, a25, a26;
0, 0, 0, 0, a35, a36;
a41, a42, a43, a44, a45, a46;
a51, a52, a53, a54, a55, a56;
0, 0, 0, 0, a65, a66];
B0 = [1; -1; 0; 0; 0; 0]; % Example: solve Ax=B for unit source
x0 = A0 B0;
x0 = x0(1);
% val0 =(-3/(4*pi*(1 + 3*x0)));
val0 =-1/(4*pi*(1 +c)*x0);
endplot after taking inversion Laplace to val function and then divide it on val0 function for steady case, it should give 1 as t tends to infty
function main_calculation()
clc; close all;
t_values = logspace(0, 2, 100);
beta1_vals = [0.001, 1, 3, 50]; % 4 specific curves to compare
params = struct(‘c’,0.0000001, ‘s1’,0.00001, ‘s2’,0.1, ‘rho1’,0, ‘beta2’,10000, ‘b’,2,’j’,0.1);
figure; hold on;
colors = [‘r’, ‘g’, ‘b’, ‘m’]; % Three colors for the three curves
for k = 1:length(beta1_vals)
params.beta1 = beta1_vals(k);
u_s1 = zeros(size(t_values));
u_steady = zeros(size(t_values));
for i = 1:length(t_values)
% Compute both inversions at the current beta1
u_s1(i) = talbot_inversion(@(s) U_func_s1(s, params), t_values(i));
u_steady(i) = talbot_inversion(@(s) U_func_steady(s, params), t_values(i));
end
% Calculate ratio
ratio = u_s1 ./ (u_steady + eps);
% Plot
plot(t_values, ratio, ‘Color’, colors(k), ‘LineWidth’, 2, …
‘DisplayName’, sprintf(‘\beta_1=%.2f’, beta1_vals(k)));
end
set(gca, ‘XScale’, ‘log’);
xlabel(‘Time (t)’); ylabel(‘Ratio (U_{osc}/U_{std})’);
legend(‘Location’, ‘best’); grid on;
% title(‘Effect of Micropolar Coupling on Dispersion Ratio’);
end
% — Rename your existing logic to avoid conflicts —
function val = U_func_s1(sigma, p)
% Copy your logic from the original calc_s1.m U_function here
% Ensure it returns ‘val’
% Unpack
b = p.b; c = p.c; s1 = p.s1; s2 = p.s2;
beta1 = p.beta1; beta2 = p.beta2; j = p.j; rho1 = p.rho1;
% Define S (sum) and P (product)
S = (4*c*s1 + (1 + j*s1*(c + 1))*sigma) / (c + 1);
P = (sigma*s1*(4*c + j*sigma)) / (c + 1);
% Calculate alpha1 and alpha2
alpha1 = sqrt((S – sqrt(S^2 – 4*P)) / 2);
alpha2 = sqrt((S + sqrt(S^2 – 4*P)) / 2);
%A1,B1,C1,D1,E1,F1
a11= 2; a12= 2 ;
a13= 2*besselk(0.3e1 / 0.2e1, alpha1);a14= 2 * besselk(0.3e1 / 0.2e1, alpha2) ;
a15= 2*besseli(0.3e1 / 0.2e1, alpha1);a16= 2 * besseli(0.3e1 / 0.2e1, alpha2);
a21= -(beta1 + 4 * c + 2) / beta1 ;
a22= 2 * (beta1 – 2 * c + 2) / beta1 ;
a23=(-alpha1 * (beta1 + 2) / beta1 * besselk(0.1e1 / 0.2e1, alpha1) – (beta1 + 4 * c + sigma + 2) / beta1 * besselk(0.3e1 / 0.2e1, alpha1));
a24=(-alpha2 * (beta1 + 2) / beta1 * besselk(0.1e1 / 0.2e1, alpha2) – (beta1 + 4 * c + sigma + 2) / beta1 * besselk(0.3e1 / 0.2e1, alpha2));
a25=(alpha1 * (beta1 + 2) / beta1 * besseli(0.1e1 / 0.2e1, alpha1) – (beta1 + 4 * c + sigma + 2) / beta1 * besseli(0.3e1 / 0.2e1, alpha1)) ;
a26=(alpha2 * (beta1 + 2) / beta1 * besseli(0.1e1 / 0.2e1, alpha2) – (beta1 + 4 * c + sigma + 2) / beta1 * besseli(0.3e1 / 0.2e1, alpha2));
%a31=0;a32=0;
a33= ((alpha1 ^ 2 * beta2 * c * s1 * s2 + alpha1 ^ 2 * beta2 * s1 * s2 – alpha1 ^ 2 * c * s1 – beta2 * s1 * s2 * sigma + alpha1 ^ 2 * c – alpha1 ^ 2 * s1 + alpha1 ^ 2 + s1 * sigma – sigma) / beta2 / s1 / c / s2 * besselk(0.3e1 / 0.2e1, alpha1) / 0.2e1 + 0.1e1 / beta2 / s2 / c * (alpha1 ^ 2 * c + alpha1 ^ 2 – sigma) * besselk(0.5e1 / 0.2e1, alpha1) * alpha1 / 0.2e1) ;
a34= ((alpha2 ^ 2 * beta2 * c * s1 * s2 + alpha2 ^ 2 * beta2 * s1 * s2 – alpha2 ^ 2 * c * s1 – beta2 * s1 * s2 * sigma + alpha2 ^ 2 * c – alpha2 ^ 2 * s1 + alpha2 ^ 2 + s1 * sigma – sigma) / beta2 / s1 / c / s2 * besselk(0.3e1 / 0.2e1, alpha2) / 0.2e1 + 0.1e1 / beta2 / s2 / c * (alpha2 ^ 2 * c + alpha2 ^ 2 – sigma) * besselk(0.5e1 / 0.2e1, alpha2) * alpha2 / 0.2e1) ;
a35= ((alpha1 ^ 2 * beta2 * c * s1 * s2 + alpha1 ^ 2 * beta2 * s1 * s2 – alpha1 ^ 2 * c * s1 – beta2 * s1 * s2 * sigma + alpha1 ^ 2 * c – alpha1 ^ 2 * s1 + alpha1 ^ 2 + s1 * sigma – sigma) / beta2 / s1 / c / s2 * besseli(0.3e1 / 0.2e1, alpha1) / 0.2e1 – 0.1e1 / beta2 / s2 / c * (alpha1 ^ 2 * c + alpha1 ^ 2 – sigma) * alpha1 * besseli(0.5e1 / 0.2e1, alpha1) / 0.2e1) ;
a36= ((alpha2 ^ 2 * beta2 * c * s1 * s2 + alpha2 ^ 2 * beta2 * s1 * s2 – alpha2 ^ 2 * c * s1 – beta2 * s1 * s2 * sigma + alpha2 ^ 2 * c – alpha2 ^ 2 * s1 + alpha2 ^ 2 + s1 * sigma – sigma) / beta2 / s1 / c / s2 * besseli(0.3e1 / 0.2e1, alpha2) / 0.2e1 – 0.1e1 / beta2 / s2 / c * (alpha2 ^ 2 * c + alpha2 ^ 2 – sigma) * besseli(0.5e1 / 0.2e1, alpha2) * alpha2 / 0.2e1) ;
a41=2 / b ^ 3; a42= 2;
a43= 2 * b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha1) ;
a44= 2 * b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha2);
a45= 2 * b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha1) ;
a46= 2 * b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha2);
a51= -1 / b ^ 3 ;a52= 2 ;
a53= (-b ^ (-0.1e1 / 0.2e1) * alpha1 * besselk(0.1e1 / 0.2e1, b * alpha1) – b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha1)) ;
a54= (-b ^ (-0.1e1 / 0.2e1) * alpha2 * besselk(0.1e1 / 0.2e1, b * alpha2) – b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha2)) ;
a55= (b ^ (-0.1e1 / 0.2e1) * alpha1 * besseli(0.1e1 / 0.2e1, b * alpha1) – b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha1)) ;
a56= (b ^ (-0.1e1 / 0.2e1) * alpha2 * besseli(0.1e1 / 0.2e1, b * alpha2) – b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha2)) ;
%a61=0; a62=0;
a63= 0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha1 ^ 2 * b ^ 2 * c + alpha1 ^ 2 * b ^ 2 – b ^ 2 * sigma) * besselk(0.3e1 / 0.2e1, b * alpha1)/ 0.2e1;
a64= 0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha2 ^ 2 * b ^ 2 * c + alpha2 ^ 2 * b ^ 2 – b ^ 2 * sigma) * besselk(0.3e1 / 0.2e1, b * alpha2) / 0.2e1;
a65= 0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha1 ^ 2 * b ^ 2 * c + alpha1 ^ 2 * b ^ 2 – b ^ 2 * sigma) * besseli(0.3e1 / 0.2e1, b * alpha1) / 0.2e1 ;
a66= 0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha2 ^ 2 * b ^ 2 * c + alpha2 ^ 2 * b ^ 2 – b ^ 2 * sigma) * besseli(0.3e1 / 0.2e1, b * alpha2) / 0.2e1;
% Construct the 6×6 matrix A and column vector B
A = [a11, a12, a13, a14, a15, a16;
a21, a22, a23, a24, a25, a26;
0, 0, a33, a34, a35, a36;
a41, a42, a43, a44, a45, a46;
a51, a52, a53, a54, a55, a56;
0, 0, a63, a64, a65, a66];
B = [1; -1; 0; 0; 0; 0]; % Example: solve Ax=B for unit source
x = A B;
x1 = x(1);
val =( 3/(4*pi*sigma^2*(rho1 – 1 – 3*x1)));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function val0 = U_func_steady(sigma, p)
% Copy your logic from the original calc_steady.m U_function here
% Ensure it returns ‘val0’
% Unpack
b = p.b; c = p.c; s1 = p.s1; s2 = p.s2;
beta1 = p.beta1; beta2 = p.beta2;
% Define S (sum) and P (product)
alpha = sqrt(4*c*s1 /(c + 1));
a11= 2; a12= 2 ; a13= 2 ; a14=2; a15=2 * besselk(0.3e1 / 0.2e1, alpha); a16=2 * besseli(0.3e1 / 0.2e1, alpha);
a21 = -(beta1 + 4 * c + 2) / beta1 ;
a22=2 * (beta1 – 2 * c + 2) / beta1 ;
a23=2 * (2 * beta1 – 7 * c – 1) / beta1;
a24=(beta1 – 2 * c + 4) / beta1 ;
a25=(-alpha * (beta1 + 2) / beta1 * besselk(0.1e1 / 0.2e1, alpha) – (beta1 + 4 * c + sigma + 2) / beta1 * besselk(0.3e1 / 0.2e1, alpha));
a26=(alpha * (beta1 + 2) / beta1 * besseli(0.1e1 / 0.2e1, alpha) – (beta1 + 4 * c + sigma + 2) / beta1 * besseli(0.3e1 / 0.2e1, alpha));
%a31=0;a32=0;a33=0;a34=0;
a35= ((alpha ^ 2 * beta2 * c * s1 * s2 + alpha ^ 2 * beta2 * s1 * s2 – alpha ^ 2 * c * s1 – beta2 * s1 * s2 * sigma + alpha ^ 2 * c – alpha ^ 2 * s1 + alpha ^ 2 + s1 * sigma – sigma) / beta2 / s1 / c / s2 * besselk(0.3e1 / 0.2e1, alpha) / 0.2e1 + 0.1e1 / beta2 / s2 / c * (alpha ^ 2 * c + alpha ^ 2 – sigma) * besselk(0.5e1 / 0.2e1, alpha) * alpha / 0.2e1) ;
a36= ((alpha ^ 2 * beta2 * c * s1 * s2 + alpha ^ 2 * beta2 * s1 * s2 – alpha ^ 2 * c * s1 – beta2 * s1 * s2 * sigma + alpha ^ 2 * c – alpha ^ 2 * s1 + alpha ^ 2 + s1 * sigma – sigma) / beta2 / s1 / c / s2 * besseli(0.3e1 / 0.2e1, alpha) / 0.2e1 – 0.1e1 / beta2 / s2 / c * (alpha ^ 2 * c + alpha ^ 2 – sigma) * besseli(0.5e1 / 0.2e1, alpha) * alpha / 0.2e1) ;
a41=2 / b ^ 3; a42=2 ; a43=2 * b ^ 2 ; a44=2 / b ;
a45=2 * b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha);
a46=2 * b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha);
a51= -1 / b ^ 3; a52= 2 ; a53=4 * b ^ 2 ;a54=1 / b ;
a55=(-b ^ (-0.1e1 / 0.2e1) * alpha * besselk(0.1e1 / 0.2e1, b * alpha) – b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha)) ;
a56=(b ^ (-0.1e1 / 0.2e1) * alpha * besseli(0.1e1 / 0.2e1, b * alpha) – b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha)) ;
%a61=0;a62=0;a63=0;a64=0;
a65= 0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha ^ 2 * b ^ 2 * c + alpha ^ 2 * b ^ 2 – b ^ 2 * sigma) * besselk(0.3e1 / 0.2e1, b * alpha) / 0.2e1;
a66= 0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha ^ 2 * b ^ 2 * c + alpha ^ 2 * b ^ 2 – b ^ 2 * sigma) * besseli(0.3e1 / 0.2e1, b * alpha) / 0.2e1;
% Construct the 6×6 matrix A and column vector B
A0 = [a11, a12, a13, a14, a15, a16;
a21, a22, a23, a24, a25, a26;
0, 0, 0, 0, a35, a36;
a41, a42, a43, a44, a45, a46;
a51, a52, a53, a54, a55, a56;
0, 0, 0, 0, a65, a66];
B0 = [1; -1; 0; 0; 0; 0]; % Example: solve Ax=B for unit source
x0 = A0 B0;
x0 = x0(1);
% val0 =(-3/(4*pi*(1 + 3*x0)));
val0 =-1/(4*pi*(1 +c)*x0);
end plot after taking inversion Laplace to val function and then divide it on val0 function for steady case, it should give 1 as t tends to infty
function main_calculation()
clc; close all;
t_values = logspace(0, 2, 100);
beta1_vals = [0.001, 1, 3, 50]; % 4 specific curves to compare
params = struct(‘c’,0.0000001, ‘s1’,0.00001, ‘s2’,0.1, ‘rho1’,0, ‘beta2’,10000, ‘b’,2,’j’,0.1);
figure; hold on;
colors = [‘r’, ‘g’, ‘b’, ‘m’]; % Three colors for the three curves
for k = 1:length(beta1_vals)
params.beta1 = beta1_vals(k);
u_s1 = zeros(size(t_values));
u_steady = zeros(size(t_values));
for i = 1:length(t_values)
% Compute both inversions at the current beta1
u_s1(i) = talbot_inversion(@(s) U_func_s1(s, params), t_values(i));
u_steady(i) = talbot_inversion(@(s) U_func_steady(s, params), t_values(i));
end
% Calculate ratio
ratio = u_s1 ./ (u_steady + eps);
% Plot
plot(t_values, ratio, ‘Color’, colors(k), ‘LineWidth’, 2, …
‘DisplayName’, sprintf(‘\beta_1=%.2f’, beta1_vals(k)));
end
set(gca, ‘XScale’, ‘log’);
xlabel(‘Time (t)’); ylabel(‘Ratio (U_{osc}/U_{std})’);
legend(‘Location’, ‘best’); grid on;
% title(‘Effect of Micropolar Coupling on Dispersion Ratio’);
end
% — Rename your existing logic to avoid conflicts —
function val = U_func_s1(sigma, p)
% Copy your logic from the original calc_s1.m U_function here
% Ensure it returns ‘val’
% Unpack
b = p.b; c = p.c; s1 = p.s1; s2 = p.s2;
beta1 = p.beta1; beta2 = p.beta2; j = p.j; rho1 = p.rho1;
% Define S (sum) and P (product)
S = (4*c*s1 + (1 + j*s1*(c + 1))*sigma) / (c + 1);
P = (sigma*s1*(4*c + j*sigma)) / (c + 1);
% Calculate alpha1 and alpha2
alpha1 = sqrt((S – sqrt(S^2 – 4*P)) / 2);
alpha2 = sqrt((S + sqrt(S^2 – 4*P)) / 2);
%A1,B1,C1,D1,E1,F1
a11= 2; a12= 2 ;
a13= 2*besselk(0.3e1 / 0.2e1, alpha1);a14= 2 * besselk(0.3e1 / 0.2e1, alpha2) ;
a15= 2*besseli(0.3e1 / 0.2e1, alpha1);a16= 2 * besseli(0.3e1 / 0.2e1, alpha2);
a21= -(beta1 + 4 * c + 2) / beta1 ;
a22= 2 * (beta1 – 2 * c + 2) / beta1 ;
a23=(-alpha1 * (beta1 + 2) / beta1 * besselk(0.1e1 / 0.2e1, alpha1) – (beta1 + 4 * c + sigma + 2) / beta1 * besselk(0.3e1 / 0.2e1, alpha1));
a24=(-alpha2 * (beta1 + 2) / beta1 * besselk(0.1e1 / 0.2e1, alpha2) – (beta1 + 4 * c + sigma + 2) / beta1 * besselk(0.3e1 / 0.2e1, alpha2));
a25=(alpha1 * (beta1 + 2) / beta1 * besseli(0.1e1 / 0.2e1, alpha1) – (beta1 + 4 * c + sigma + 2) / beta1 * besseli(0.3e1 / 0.2e1, alpha1)) ;
a26=(alpha2 * (beta1 + 2) / beta1 * besseli(0.1e1 / 0.2e1, alpha2) – (beta1 + 4 * c + sigma + 2) / beta1 * besseli(0.3e1 / 0.2e1, alpha2));
%a31=0;a32=0;
a33= ((alpha1 ^ 2 * beta2 * c * s1 * s2 + alpha1 ^ 2 * beta2 * s1 * s2 – alpha1 ^ 2 * c * s1 – beta2 * s1 * s2 * sigma + alpha1 ^ 2 * c – alpha1 ^ 2 * s1 + alpha1 ^ 2 + s1 * sigma – sigma) / beta2 / s1 / c / s2 * besselk(0.3e1 / 0.2e1, alpha1) / 0.2e1 + 0.1e1 / beta2 / s2 / c * (alpha1 ^ 2 * c + alpha1 ^ 2 – sigma) * besselk(0.5e1 / 0.2e1, alpha1) * alpha1 / 0.2e1) ;
a34= ((alpha2 ^ 2 * beta2 * c * s1 * s2 + alpha2 ^ 2 * beta2 * s1 * s2 – alpha2 ^ 2 * c * s1 – beta2 * s1 * s2 * sigma + alpha2 ^ 2 * c – alpha2 ^ 2 * s1 + alpha2 ^ 2 + s1 * sigma – sigma) / beta2 / s1 / c / s2 * besselk(0.3e1 / 0.2e1, alpha2) / 0.2e1 + 0.1e1 / beta2 / s2 / c * (alpha2 ^ 2 * c + alpha2 ^ 2 – sigma) * besselk(0.5e1 / 0.2e1, alpha2) * alpha2 / 0.2e1) ;
a35= ((alpha1 ^ 2 * beta2 * c * s1 * s2 + alpha1 ^ 2 * beta2 * s1 * s2 – alpha1 ^ 2 * c * s1 – beta2 * s1 * s2 * sigma + alpha1 ^ 2 * c – alpha1 ^ 2 * s1 + alpha1 ^ 2 + s1 * sigma – sigma) / beta2 / s1 / c / s2 * besseli(0.3e1 / 0.2e1, alpha1) / 0.2e1 – 0.1e1 / beta2 / s2 / c * (alpha1 ^ 2 * c + alpha1 ^ 2 – sigma) * alpha1 * besseli(0.5e1 / 0.2e1, alpha1) / 0.2e1) ;
a36= ((alpha2 ^ 2 * beta2 * c * s1 * s2 + alpha2 ^ 2 * beta2 * s1 * s2 – alpha2 ^ 2 * c * s1 – beta2 * s1 * s2 * sigma + alpha2 ^ 2 * c – alpha2 ^ 2 * s1 + alpha2 ^ 2 + s1 * sigma – sigma) / beta2 / s1 / c / s2 * besseli(0.3e1 / 0.2e1, alpha2) / 0.2e1 – 0.1e1 / beta2 / s2 / c * (alpha2 ^ 2 * c + alpha2 ^ 2 – sigma) * besseli(0.5e1 / 0.2e1, alpha2) * alpha2 / 0.2e1) ;
a41=2 / b ^ 3; a42= 2;
a43= 2 * b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha1) ;
a44= 2 * b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha2);
a45= 2 * b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha1) ;
a46= 2 * b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha2);
a51= -1 / b ^ 3 ;a52= 2 ;
a53= (-b ^ (-0.1e1 / 0.2e1) * alpha1 * besselk(0.1e1 / 0.2e1, b * alpha1) – b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha1)) ;
a54= (-b ^ (-0.1e1 / 0.2e1) * alpha2 * besselk(0.1e1 / 0.2e1, b * alpha2) – b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha2)) ;
a55= (b ^ (-0.1e1 / 0.2e1) * alpha1 * besseli(0.1e1 / 0.2e1, b * alpha1) – b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha1)) ;
a56= (b ^ (-0.1e1 / 0.2e1) * alpha2 * besseli(0.1e1 / 0.2e1, b * alpha2) – b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha2)) ;
%a61=0; a62=0;
a63= 0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha1 ^ 2 * b ^ 2 * c + alpha1 ^ 2 * b ^ 2 – b ^ 2 * sigma) * besselk(0.3e1 / 0.2e1, b * alpha1)/ 0.2e1;
a64= 0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha2 ^ 2 * b ^ 2 * c + alpha2 ^ 2 * b ^ 2 – b ^ 2 * sigma) * besselk(0.3e1 / 0.2e1, b * alpha2) / 0.2e1;
a65= 0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha1 ^ 2 * b ^ 2 * c + alpha1 ^ 2 * b ^ 2 – b ^ 2 * sigma) * besseli(0.3e1 / 0.2e1, b * alpha1) / 0.2e1 ;
a66= 0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha2 ^ 2 * b ^ 2 * c + alpha2 ^ 2 * b ^ 2 – b ^ 2 * sigma) * besseli(0.3e1 / 0.2e1, b * alpha2) / 0.2e1;
% Construct the 6×6 matrix A and column vector B
A = [a11, a12, a13, a14, a15, a16;
a21, a22, a23, a24, a25, a26;
0, 0, a33, a34, a35, a36;
a41, a42, a43, a44, a45, a46;
a51, a52, a53, a54, a55, a56;
0, 0, a63, a64, a65, a66];
B = [1; -1; 0; 0; 0; 0]; % Example: solve Ax=B for unit source
x = A B;
x1 = x(1);
val =( 3/(4*pi*sigma^2*(rho1 – 1 – 3*x1)));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function val0 = U_func_steady(sigma, p)
% Copy your logic from the original calc_steady.m U_function here
% Ensure it returns ‘val0’
% Unpack
b = p.b; c = p.c; s1 = p.s1; s2 = p.s2;
beta1 = p.beta1; beta2 = p.beta2;
% Define S (sum) and P (product)
alpha = sqrt(4*c*s1 /(c + 1));
a11= 2; a12= 2 ; a13= 2 ; a14=2; a15=2 * besselk(0.3e1 / 0.2e1, alpha); a16=2 * besseli(0.3e1 / 0.2e1, alpha);
a21 = -(beta1 + 4 * c + 2) / beta1 ;
a22=2 * (beta1 – 2 * c + 2) / beta1 ;
a23=2 * (2 * beta1 – 7 * c – 1) / beta1;
a24=(beta1 – 2 * c + 4) / beta1 ;
a25=(-alpha * (beta1 + 2) / beta1 * besselk(0.1e1 / 0.2e1, alpha) – (beta1 + 4 * c + sigma + 2) / beta1 * besselk(0.3e1 / 0.2e1, alpha));
a26=(alpha * (beta1 + 2) / beta1 * besseli(0.1e1 / 0.2e1, alpha) – (beta1 + 4 * c + sigma + 2) / beta1 * besseli(0.3e1 / 0.2e1, alpha));
%a31=0;a32=0;a33=0;a34=0;
a35= ((alpha ^ 2 * beta2 * c * s1 * s2 + alpha ^ 2 * beta2 * s1 * s2 – alpha ^ 2 * c * s1 – beta2 * s1 * s2 * sigma + alpha ^ 2 * c – alpha ^ 2 * s1 + alpha ^ 2 + s1 * sigma – sigma) / beta2 / s1 / c / s2 * besselk(0.3e1 / 0.2e1, alpha) / 0.2e1 + 0.1e1 / beta2 / s2 / c * (alpha ^ 2 * c + alpha ^ 2 – sigma) * besselk(0.5e1 / 0.2e1, alpha) * alpha / 0.2e1) ;
a36= ((alpha ^ 2 * beta2 * c * s1 * s2 + alpha ^ 2 * beta2 * s1 * s2 – alpha ^ 2 * c * s1 – beta2 * s1 * s2 * sigma + alpha ^ 2 * c – alpha ^ 2 * s1 + alpha ^ 2 + s1 * sigma – sigma) / beta2 / s1 / c / s2 * besseli(0.3e1 / 0.2e1, alpha) / 0.2e1 – 0.1e1 / beta2 / s2 / c * (alpha ^ 2 * c + alpha ^ 2 – sigma) * besseli(0.5e1 / 0.2e1, alpha) * alpha / 0.2e1) ;
a41=2 / b ^ 3; a42=2 ; a43=2 * b ^ 2 ; a44=2 / b ;
a45=2 * b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha);
a46=2 * b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha);
a51= -1 / b ^ 3; a52= 2 ; a53=4 * b ^ 2 ;a54=1 / b ;
a55=(-b ^ (-0.1e1 / 0.2e1) * alpha * besselk(0.1e1 / 0.2e1, b * alpha) – b ^ (-0.3e1 / 0.2e1) * besselk(0.3e1 / 0.2e1, b * alpha)) ;
a56=(b ^ (-0.1e1 / 0.2e1) * alpha * besseli(0.1e1 / 0.2e1, b * alpha) – b ^ (-0.3e1 / 0.2e1) * besseli(0.3e1 / 0.2e1, b * alpha)) ;
%a61=0;a62=0;a63=0;a64=0;
a65= 0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha ^ 2 * b ^ 2 * c + alpha ^ 2 * b ^ 2 – b ^ 2 * sigma) * besselk(0.3e1 / 0.2e1, b * alpha) / 0.2e1;
a66= 0.1e1 / c * b ^ (-0.5e1 / 0.2e1) * (alpha ^ 2 * b ^ 2 * c + alpha ^ 2 * b ^ 2 – b ^ 2 * sigma) * besseli(0.3e1 / 0.2e1, b * alpha) / 0.2e1;
% Construct the 6×6 matrix A and column vector B
A0 = [a11, a12, a13, a14, a15, a16;
a21, a22, a23, a24, a25, a26;
0, 0, 0, 0, a35, a36;
a41, a42, a43, a44, a45, a46;
a51, a52, a53, a54, a55, a56;
0, 0, 0, 0, a65, a66];
B0 = [1; -1; 0; 0; 0; 0]; % Example: solve Ax=B for unit source
x0 = A0 B0;
x0 = x0(1);
% val0 =(-3/(4*pi*(1 + 3*x0)));
val0 =-1/(4*pi*(1 +c)*x0);
end plot, inversion laplce MATLAB Answers — New Questions
