0

My calculation involves cosh(x) and sinh(x) when x is around 700 - 1000 which reaches MATLAB's limit and the result is NaN. The problem in the code is elastic_restor_coeff rises when radius is small (below 5e-9 in the code). My goal is to do another integral over a radius distribution from 1e-9 to 100e-9 which is still a work in progress because I get stuck at this problem.

My work around solution right now is to approximate the real part of chi_para with a step function when threshold2 hits a value of about 300. The number 300 is obtained from using the lowest possible value of radius and look at the cut-off value from the plot. I think this approach is not good enough for actual calculation since this value changes with radius so I am looking for a better approximation method. Also, the imaginary part of chi_para is difficult to approximate since it looks like a pulse instead of a step.

Here is my code without an integration over a radius distribution.

k_B = 1.38e-23;
T = 296;
radius = [5e-9,10e-9, 20e-9, 30e-9,100e-9];
fric_coeff = 8*pi*1e-3.*radius.^3;
elastic_restor_coeff = 8*pi*1.*radius.^3;
time_const = fric_coeff/elastic_restor_coeff;
omega_ar = logspace(-6,6,60);
chi_para = zeros(1,length(omega_ar));
chi_perpen = zeros(1,length(omega_ar));
threshold = zeros(1,length(omega_ar));
threshold2 = zeros(1,length(omega_ar));
for i = 1:length(radius)
    for k = 1:length(omega_ar)
        omega = omega_ar(k);
        fric_coeff = 8*pi*1e-3.*radius(i).^3;
        elastic_restor_coeff = 8*pi*1.*radius(i).^3;
        time_const = fric_coeff/elastic_restor_coeff;
        G_para_func = @(t) ((cosh(2*k_B*T./elastic_restor_coeff.*exp(-t./time_const))-1).*exp(1i.*omega.*t))./(cosh(2*k_B*T./elastic_restor_coeff)-1);
        G_perpen_func = @(t) ((sinh(2*k_B*T./elastic_restor_coeff.*exp(-t./time_const))).*exp(1i.*omega.*t))./(sinh(2*k_B*T./elastic_restor_coeff));

        chi_para(k) = (1 + 1i*omega*integral(G_para_func, 0, inf));
        chi_perpen(k) = (1 + 1i*omega*integral(G_perpen_func, 0, inf));
        threshold(k) = 2*k_B*T./elastic_restor_coeff*omega;
        threshold2(k) = 2*k_B*T./elastic_restor_coeff*(omega*time_const - 1);
    end
    figure(1);
    semilogx(omega_ar,real(chi_para),omega_ar,imag(chi_para));
    hold on;
    figure(2);
    semilogx(omega_ar,real(chi_perpen),omega_ar,imag(chi_perpen));
    hold on;
end

Here is the simplified function that I would like to approximate:

cosh equation

where x is iterated in a loop and the maximum value of x is about 700.

Nabs
  • 362
  • 1
  • 3
  • 10
  • 1
    Have you considered simply performing the calculations involving `sinh()` and `cosh()` in logarithmic space, then exponentiate at the end? Consider the asymptotic behavior of these functions for arguments of large magnitude. It might also be worthwhile to revisit the course of computation that leads to large arguments being passed to these functions in the first place. Are alternative computational trajectories feasible? – njuffa Mar 18 '18 at 23:12
  • you mean log(cosh(x)) approximated to x + log(1/2) when x is large? That's brilliant! Never crossed my mind. However, I also have exp(-t) when t is from 0 to inf multiplied with x. So, in this case, I have to separate this integration into two parts when x*exp(-t) is large enough and when it is small that the approximation does not work. Is my thought correct? – Nabs Mar 19 '18 at 01:38
  • Yes, for |x| > 30 or so we have cosh(x) = 0.5*exp(abs(x)) to within double-precision accuracy. I don't know what exactly you are computing with the help of `sinh` and `cosh`, but it seems you are on the right track. – njuffa Mar 19 '18 at 03:06
  • Thanks a lot! It's unfortunate that stackoverflow does not support Latex so I could put up an equation. But I will try to implement this and see if it there is no other issues coming out. – Nabs Mar 19 '18 at 03:13
  • I am wondering how to actually implement this because of the integration from 0 to inf of the function. If it was calculating point by point then I could see how to do it. – Nabs Mar 19 '18 at 12:24
  • Sorry, I am a CS guy, not a mathematician. In numerical codes, a standard trick to avoid floating-point overflow and underflow in intermediate computation is to temporarily move to logarithmic space. That's all I meant to point out. – njuffa Mar 19 '18 at 15:56
  • Thank you for your suggestion anyway but I am still trying to find a way to implement this in the code. – Nabs Mar 20 '18 at 00:16

0 Answers0