# Residue failed to reproduce this residue in a paper

This problem is a more clearly state of my original problem which is off topic, so I reconsider my problem carefully based on previous discussions, and gives the following valid code to reproduce my problem with all the necessary information.

I use Mathematica to confirm a result from a paper Unsteady unidirectional flow of Bingham fluid between parallel plates with different given volume flow rate conditions which is to calculate the residue of a expression at zero.

1. Reproduce of my problem:

The following code shows my efforts, but failed to reproduce the analytical solution (resAna in the code) in the papaer.

ClearAll["Global*"]

varCons = {ν > 0, up > 0, h > h0 >= 0, h > y >= 0};
m = Sqrt[s]/Sqrt[ν];
Δ = Sinh[m h] Sinh[m h0] - Cosh[m h] Cosh[m h0];
Ξ =  Cosh[m h0] (Sinh[m h] - Sinh[m h0]) - Sinh[m h0] (Cosh[m h] - Cosh[m h0]);
Ω = (m h (Δ + Cosh[m h0] Cosh[m y] - Sinh[m h0] Sinh[m y]))/(m h Δ +
m h0 (Cosh[m h0]^2 - Sinh[m h0]^2) + Ξ);
u = FullSimplify[up/(t0 s^2) Ω E^(s t), varCons];

Residue[u, {s, 0}]
(* this failed gives the desired 'resAna' as below *)

X1 = (1/6 (h^3 - y^3) h0 + 1/6 (h - y) h0^3 - 1/4 (h^2 - y^2) h0^2 -1/24 (h^4 - y^4)) (1/ν)^(5/2);
X2 = ((h - y) h0 - 1/2 (h^2 - y^2)) t (1/ν)^(3/2);
X3 = (h^5/30 - (h^4 h0)/8 + (h^3 h0^2)/6 - (h^2 h0^3)/4 + (21 h0^5)/120) (1/2 (h^2 - y^2) - (h - y) h0)/(h0^3/6 - (h0 h^2)/2 + h^3/3) (1/ν)^(5/2);
resAna = -((up h)/t0) (X1 + X2 + X3)/((h0^3/6 - (h0 h^2)/2 + h^3/3) (1/ν)^(3/2));
(* resAna is the analytical solution in the paper *)


2. Some numerical tries to get a clue or two

NResidue seems works after given a set of numeric values, but Residue still doesn't work

t0 = 1/10;
t = 2 t0;
up = 1/100;
ν = 2/3580;
h = 1/1000;
h0 = h/8;
y = h/9;

Needs["NumericalCalculus"]
NResidue[u, {s, 0}, WorkingPrecision -> 20, PrecisionGoal -> 20] // Chop

N[resAna, 20]

Residue[u, {s, 0}]


3. I had thought that Residue should work as in this example by @xzczd

ClearAll["Global*"]
r = 1;
f[p_, ξ_] = -(5 p Sqrt[(5 p^2)/6 + ξ^2])/(4 (-4 ξ^2 Sqrt[(5 p^2)/6 + ξ^2] Sqrt[(5 p^2)/2 + ξ^2] + ((5 p^2)/2 + 2 ξ^2)^2));
poles = p /. Simplify[Solve[Denominator[f[p, ξ]] == 0, p], ξ > 0]
Residue[f[p, ξ], {p, #}] & /@ poles;
Simplify[%, ξ > 0] // N

• The last example is redundant. It's taken from this post of mine, right? The relevance between these 2 problems isn't obvious. BTW that problem isn't simple at all, I think. – xzczd Jun 17 '16 at 8:52
• @xzczd, Yes, the example that Residue works is from your post. Of course, the problem in your post isn't simple. By "simple" I meant that the syntax of Residue is simple, just to reassure that I am using Residue correctly. Sorry for the misunderstanding caused by my poor english. – xinxin guo Jun 17 '16 at 9:21

I. A summary for the failing trial

1. Residue can't calculate the residue of u directly. (It's hard to tell whether it's a bug or not, Residue never promises he will calculate every known residue, anyway.)

2. SeriesCoefficient won't give desired answer in the following case:

SeriesCoefficient[u, {s, 0, -1}]


If it gave the desired answer, we would be able to calculate the residue by finding Laurent series expansions.

3. Limit won't give desired answer in the following case:

Limit[D[s^2 u, s], s -> 0]


If it gave the desired answer, we would be able to calculate the residue with limit formula.

II. Solution

By observing the structure of undesired output of SeriesCoefficient and u

I guess that maybe it's the $\sqrt{\nu}$ term that causes trouble, so I tried adding the known constraint $\nu>0$ to SeriesCoefficient and Limit, and they give the desired answer then:

resAnaTrue =
Simplify[SeriesCoefficient[u, {s, 0, -1}, Assumptions -> ν > 0], ν > 0];

resAnaTrue2 = Limit[D[s^2 u, s], s -> 0, Assumptions -> ν > 0]
(* -(h up (h - y) (h - 2 h0 + y) (2 h^3 - 6 h^2 h0 -
2 h (2 h0^2 - 10 h0 y + 5 y^2 + 60 t ν) -
h0 (7 h0^2 - 10 h0 y + 5 y^2 + 60 t ν)))/(20 (h - h0)^2 (2 h + h0)^2 t0 ν) *)

resAnaTrue == resAnaTrue2 // Simplify
(* True *)


III. The formula given in the paper is incorrect

Let's first check the result with the parameters given by you:

para = {t0 -> 1/10, t -> 2 t0, up -> 1/100, ν -> 2/3580, h -> 1/1000, h0 -> h/8, y -> h/9};
ref = NResidue[u //. para, {s, 0}, WorkingPrecision -> 64];
rst1 = N[resAna //. para, 64];
rst2 = N[resAnaTrue //. para, 64];
ref - rst1
(* -2.9855303465924184799970408331399786805478493918575833694822*10^-7 *)
ref - rst2
(* 0.*10^-66 *)


As one can see, the precision of the new derived formula is desired, while there's a significant difference (approximately 3*10^-7) between the formula given in the paper and numeric result.

With the following set of parameter, the difference will be even clearer:

para2 = {ν -> 1/2, h0 -> 1, h -> 3, t0 -> 1/5, up -> 1/5, y -> 1/3, t -> 5};
ref = NResidue[u /. para2, {s, 0}];
rst1 = N[resAna /. para2];
rst2 = N[resAnaTrue /. para2];
ref - rst1
(* -0.653061 + 1.22133*10^-13 I *)
ref - rst2
(* 9.76996*10^-15 + 1.22133*10^-13 I *)
`

Apparently the formula in the paper is wrong.

• I really appreciate your help. :) – xinxin guo Jun 18 '16 at 11:04
• @pengfei_guo Er… just to confirm, do you know that with more than 15 reputation you can now upvote by clicking the gray triangles? – xzczd Jun 18 '16 at 15:58
• @xinxinguo I know you are a Chinese in Xi'an city via your profile. BTW, which university are you in? – xyz Jun 22 '16 at 7:29
• I'm a young lecture at Xi'an university of architecture and technology 西安建筑科技大学. – xinxin guo Jun 22 '16 at 7:38
• @xinxinguo Thanks, got it:) That's mean you are a university teacher. – xyz Jun 30 '16 at 1:49