How to solve an elliptical integral with the unknown in the integrand?

Hello everyone,

I am trying to implement the exact solution to the Jeffery-Hamel flow problem according to White (1974) [Link]. The two equations of interest are here, with and α given.

1. Solution:
2. Boundary condition:

So first, I would like to determine C by using equation (2) using the following code:

Re = 100;
alpha = 15/360*2*pi;

bc = @(f, C) 1./(((1-f).*(2/3.*Re.*alpha.*(f.*f+f)+4.*alpha.*alpha.*f + C)).^(0.5));
bcInt = @(C) integral(@(f)  bc(f, C), 0, rand(1));
bcEq = @(C) bcInt(C) - 1;
C_sol = fsolve(bcEq, rand)

However, I get the following error message:

No solution found.

fsolve stopped because the relative size of the current step is less than the
value of the step size tolerance squared, but the vector of function values
is not near zero as measured by the value of the function tolerance.

<stopping criteria details>

Does anybody have an idea what I do wrong, or if I have to take a completely different path to solve this?

Thanks a lot!

Answer 1

[Alan Stevens]

Here's how to use fzero to find C. However, with Re = 100 and alpha = 15deg, the boundary condition integral never reaches 1 for any C. Reduce Re to 20, say, or alpha to 1.5 and you get a result.

C0 = 5; % initial guess
C = fzero(@bc ,C0);
disp(C)

function Z = bc(C)
    Re = 20;          % The bc never reaches 1 with Re = 100 and alpha = 15 degrees
    alpha = deg2rad(15);   
    fn = @(f) 1./(  ((1-f).*(2/3.*Re.*alpha.*(f.*f+f)+4.*alpha.*alpha.*f + C)).^(0.5)  );
    Z = 1-integral(fn,0,1;     
end

Answer 2

[David Hill]

If you plot for C=0:100, you find that bcInt never exceeds 1.

Re = 100;
alpha = 15/360*2*pi;
x=0:100;
bc = @(f, C) 1./(((1-f).*(2/3.*Re.*alpha.*(f.*f+f)+4.*alpha.*alpha.*f + C)).^(0.5));
bcInt = arrayfun(@(C) integral(@(f)  bc(f, C), 0, 1),x);%not sure why you are changing your integration limit (rand(1))
plot(x,bcInt);

Therefore, there will be no solution when you subtract 1, but if you subtract something in the range

Re = 100;
alpha = 15/360*2*pi;
bc = @(f, C) 1./(((1-f).*(2/3.*Re.*alpha.*(f.*f+f)+4.*alpha.*alpha.*f + C)).^(0.5));
bcInt = @(x)arrayfun(@(C) integral(@(f)  bc(f, C), 0, 1)-.29,x);%subtract .29
C_sol = fzero(bcInt, 100);%finds the C @.29