How to avoid inf/inf numerically for hyperbolic functions

I have an expression which consist of hyperbolic functions so to avoid the numerically instability how i can write them in exponential form in coding. I have those expressions in fortan code. Now i want to use them in MATLAB for my problem but i did not get that required form. The parameters i have are h=200, d= 38 , \gamma is the positive roots for bessel functions.

and

Now the fortan code they write for C_n is:

kmR(iM)=BesJzero(iM)                               ! this is the \gamma they just used the scaled one but analytically this \gamma.
    kmH(iM)=kmR(iM)*HsR                            ! this is \gamm*h 
    exp1kmD(iM)=exp(-kmH(iM)*dsD/wtH)              ! this is exp(-\gamma*d)
    exp2kmD(iM)=exp(-(kmH(iM)+kmH(iM))*dsD/wtH)    ! this is exp(-2 \gamma *d)
    exp2kmH(iM)=exp(-(kmH(iM)+kmH(iM)))            ! this is exp(-2 \gamma *h)
    exp2kmHmD(iM)=exp(-(kmH(iM)+kmH(iM))*HmD)      ! this is exp(-2 \gamma(h-d))

    BBm(iM)=(nuR-kmR(iM))*exp2kmHmD(iM)/exp2kmD(iM)   ! this nuR is f^2 and 
    BBm(iM)=BBm(iM)+(nuR+kmR(iM))*(1.0D0+2.0D0*exp2kmHmD(iM))
    BBm(iM)=BBm(iM)/(nuR-kmR(iM)+(nuR+kmR(iM))*exp2kmD(iM))
    BBm(iM)=BBm(iM)/(1.0D0+exp2kmHmD(iM))
    
    CCm(iM)=2.0D0*(BBm(iM)*exp2kmD(iM)-1.0D0)

And for phi_U is

CCmc(iM)=1.0D0+exp2kmHmD(iM)/exp2kmD(iM)+2.0D0*exp2kmHmD(iM)
    CCmc(iM)=CCmc(iM)/(1.0D0+exp2kmHmD(iM))
    CCmc(iM)=BBm(iM)*(1.0D0+exp2kmD(iM))-CCmc(iM)
    CCmc(iM)=CCmc(iM)*exp1kmD(iM)

Answer 1

[Walter Roberson]

Unless I have made a mistake, your expression simplifies a lot.

syms gamma__l(n) h f d
part1a = gamma__l(n)*cosh(gamma__l(n)*h)-f^2*sinh(gamma__l(n)*h);
part1b = f^2*cosh(gamma__l(n)*d)-gamma__l(n)*sinh(gamma__l(n)*d);
part2 = 1/cosh(gamma__l(n)*(h-d));
C__l(n) = part1a/part1b * part2

C__l(n) = 

part3 = C__l(n) * cosh(gamma__l(n)*d);
part4 = sinh(gamma__l(n)*h) / cosh(gamma__l(n)*(h-d));
phi(n) = part3 + part4

phi(n) = 

temp1 = simplify(phi, 'steps', 50)

temp1(n) = 

●

phi_U__tilde__l = symsum(temp1(n), n, 1, inf)

phi_U__tilde__l = 

temp2 = simplify(phi_U__tilde__l, 'steps', 50)

temp2 = 

rewrite(temp2, 'exp')

ans = 

The bottom is the expression written in terms of exponentials.

Caution: the great majority of the time, cosh() and sinh() are more accurate than writing in terms of exponentials, and have less overflow. Writing in terms of exponentials is the wrong thing to do the majority of the time.

[Javeria]

@Walter Roberson I used the hyperbolic one and it then give me NaN. So then i have this code but it is in fortan i tried alot to make the expression like they have but i couldn't.

[Walter Roberson]

The positive roots of the bessel function grow indefinitely, tending towards being πapart. You are muitiplying those by a positive constant. It doesn't take long until the product of the constant and the roots exceeds 710. At that point, cosh() of the product becomes infinite if you are using double precision.

How quickly? Well, you are multiplying by 200 in one case, and 200*pi is approaching 710, so you only get as far as the first root, since the second root is 5.52007811028631 and 5.5*200 > 710, cosh(710) --> inf.

You should therefor expect extreme numeric problems if you are using double precision. Using the exponential form instead of cosh() form will do nothing to solve the problem.

You pretty much need to switch to the Symbolic Toolbox to get anywhere.

Answer 2

[Stephen23]

1) Fortran-style exponential factoring (stable; no direct cosh/sinh)

This is a straight translation of your Fortran into MATLAB, it computes the same intermediate arrays.

h = 200;
d = 38;
N = 7;
gammaV = rand(1,3); % yourGammaVector
f = rand(1,3);

Cn = stable_coeffs_from_exp(gammaV, f, h, d)

Cn = 3×3

-2.0187 -1.9803 -1.9805 -2.0158 -1.9831 -1.9833 -2.0000 -2.0000 -2.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

function [Cn, phiU_term, BBm] = stable_coeffs_from_exp(gamma, f, h, d)
% Numerically stable coefficients using only exponentials with negative exponents.
% Inputs:
%   gamma : vector of positive Bessel zeros (γ_n^ℓ)
%   f     : scalar or vector (same size as gamma allowed); nuR = f.^2
%   h,d   : scalars
% Outputs:
%   Cn        : equals CCm in your Fortran (i.e., C_n^ℓ)
%   phiU_term : equals CCmc in your Fortran (term used in \tilde{φ}_U^ℓ)
%   BBm       : intermediate as in Fortran (often used elsewhere)

gamma = gamma(:);
nuR   = f.^2;
% Exponentials (all <= 1 for γ,h,d > 0)
exp1kmD   = exp(-gamma.*d);             % e^{-γ d}
exp2kmD   = exp(-2*gamma.*d);           % e^{-2γ d}
exp2kmH   = exp(-2*gamma.*h);           %#ok<NASGU> % kept for completeness (mirrors Fortran)
exp2kmHmD = exp(-2*gamma.*(h - d));     % e^{-2γ (h-d)}
%
% -------- Fortran's BBm --------
BBm  = (nuR - gamma).*exp2kmHmD ./ exp2kmD;
BBm  = BBm + (nuR + gamma).*(1 + 2*exp2kmHmD);
BBm  = BBm ./ (nuR - gamma + (nuR + gamma).*exp2kmD);
BBm  = BBm ./ (1 + exp2kmHmD);
%
% -------- Cn (CCm in Fortran) --------
Cn = 2*(BBm.*exp2kmD - 1);
%
% -------- phiU_term (CCmc in Fortran) --------
phiU_term  = 1 + exp2kmHmD./exp2kmD + 2*exp2kmHmD;
phiU_term  = phiU_term ./ (1 + exp2kmHmD);
phiU_term  = BBm.*(1 + exp2kmD) - phiU_term;
phiU_term  = phiU_term .* exp1kmD;
%
end

2) Direct hyperbolic form, rewritten with exponentials (also stable)

If you want to stick to the analytical formula

Cnℓ=γcosh⁡(γh)−f2sinh⁡(γh)f2cosh⁡(γd)−γsinh⁡(γd)⋅1cosh⁡(γ(h−d)),Cnℓ​=f2cosh(γd)−γsinh(γd)γcosh(γh)−f2sinh(γh)​⋅cosh(γ(h−d))1​,

do not compute exp(gamma*...) directly (it can overflow).

Compute t=e−γxt=e−γx and form cosh/sinh from tt and 1/t1/t:

Cn = Cn_from_hyperbolic(gammaV, f, h, d)

Cn = 3×3

-2.0187 -1.9803 -1.9805 -2.0158 -1.9831 -1.9833 -2.0000 -2.0000 -2.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

function Cn = Cn_from_hyperbolic(gamma, f, h, d)
% C_n^ℓ computed from the hyperbolic formula using only exp(-γx) terms.
gamma = gamma(:);
t_h  = exp(-gamma.*h);          % e^{-γ h}   (<= 1)
t_d  = exp(-gamma.*d);          % e^{-γ d}
t_hd = exp(-gamma.*(h - d));    % e^{-γ (h-d)}
% cosh(x) = (e^x + e^{-x})/2 = (1/t + t)/2  when t = e^{-x}
% sinh(x) = (e^x - e^{-x})/2 = (1/t - t)/2
cosh_h  = 0.5*(t_h.^-1 + t_h);
sinh_h  = 0.5*(t_h.^-1 - t_h);
cosh_d  = 0.5*(t_d.^-1 + t_d);
sinh_d  = 0.5*(t_d.^-1 - t_d);
cosh_hd = 0.5*(t_hd.^-1 + t_hd);
Cn = (gamma.*cosh_h - f.^2.*sinh_h ) ./ ( f.^2.*cosh_d - gamma.*sinh_d ) ./ cosh_hd;

end

Why this is stable: all exponentials are of the form exp(-γ*positive) so they’re in (0,1](0,1]. You never form exp(+γx) directly, which is what overflows for large γ*h or γ*d.

[Javeria]

@Stephen23 hello i used your direct hyperbolic form written with exponential form , then this one give me inf.

gamma=3.61969011342704;
h=200;
t_h  = exp(-gamma.*h)

t_h = 3.9601e-315

cosh_h  = 0.5*(t_h.^-1 + t_h)

cosh_h = Inf

[Torsten]

Then you should use the method from the Fortran Code. It seems the authors knew what they did.

Copied from @Stephen23 ' s code:

gamma = 3.61969011342704;
f = ...;
d = 38;
h = 200;
nuR   = f.^2;
% Exponentials (all <= 1 for γ,h,d > 0)
exp1kmD   = exp(-gamma.*d);             % e^{-γ d}
exp2kmD   = exp(-2*gamma.*d);           % e^{-2γ d}
exp2kmH   = exp(-2*gamma.*h);           %#ok<NASGU> % kept for completeness (mirrors Fortran)
exp2kmHmD = exp(-2*gamma.*(h - d));     % e^{-2γ (h-d)}
%
% -------- Fortran's BBm --------
BBm  = (nuR - gamma).*exp2kmHmD ./ exp2kmD;
BBm  = BBm + (nuR + gamma).*(1 + 2*exp2kmHmD);
BBm  = BBm ./ (nuR - gamma + (nuR + gamma).*exp2kmD);
BBm  = BBm ./ (1 + exp2kmHmD);
%
% -------- Cn (CCm in Fortran) --------
Cn = 2*(BBm.*exp2kmD - 1);
%
% -------- phiU_term (CCmc in Fortran) --------
phiU_term  = 1 + exp2kmHmD./exp2kmD + 2*exp2kmHmD;
phiU_term  = phiU_term ./ (1 + exp2kmHmD);
phiU_term  = BBm.*(1 + exp2kmD) - phiU_term;
phiU_term  = phiU_term .* exp1kmD;

[Torsten]

Here is the "proof" that both terms (Cn and phiU_term) are equal to the expressions that you want to set in your code:

syms gamma f d h
nuR   = f.^2;
% Exponentials (all <= 1 for γ,h,d > 0)
exp1kmD   = exp(-gamma.*d);             % e^{-γ d}
exp2kmD   = exp(-2*gamma.*d);           % e^{-2γ d}
exp2kmH   = exp(-2*gamma.*h);           %#ok<NASGU> % kept for completeness (mirrors Fortran)
exp2kmHmD = exp(-2*gamma.*(h - d));     % e^{-2γ (h-d)}
%
% -------- Fortran's BBm --------
BBm  = (nuR - gamma).*exp2kmHmD ./ exp2kmD;
BBm  = BBm + (nuR + gamma).*(1 + 2*exp2kmHmD);
BBm  = BBm ./ (nuR - gamma + (nuR + gamma).*exp2kmD);
BBm  = BBm ./ (1 + exp2kmHmD);
%
% -------- Cn (CCm in Fortran) --------
Cn = 2*(BBm.*exp2kmD - 1);
%
% -------- phiU_term (CCmc in Fortran) --------
phiU_term  = 1 + exp2kmHmD./exp2kmD + 2*exp2kmHmD;
phiU_term  = phiU_term ./ (1 + exp2kmHmD);
phiU_term  = BBm.*(1 + exp2kmD) - phiU_term;
phiU_term  = phiU_term .* exp1kmD;

%Check expressions for equality
cosh_gammah = (exp(gamma*h)+exp(-gamma*h))/2;
sinh_gammah = (exp(gamma*h)-exp(-gamma*h))/2;
cosh_gammad = (exp(gamma*d)+exp(-gamma*d))/2;
sinh_gammad = (exp(gamma*d)-exp(-gamma*d))/2;
cosh_gammahmd = (exp(gamma*(h-d))+exp(-gamma*(h-d)))/2;
Cn1 = (gamma*cosh_gammah-nuR*sinh_gammah)/(nuR*cosh_gammad-gamma*sinh_gammad)/cosh_gammahmd;
phiU_term1 = Cn1*cosh_gammad+sinh_gammah/cosh_gammahmd;
simplify(Cn-Cn1,'steps',50)

ans = 

0

simplify(phiU_term-phiU_term1,'steps',50)

ans = 

0