Help using fzero to solve a specific equation

Hello everyone,

I'm in dire need to solve the following equation, the von Karman equation:

1/sqrt(f) = 4*log10(Re*sqrt(f)) - 0.4

Re is given by:

Min_Re = input('Choose minimum value of Re: ');
Max_Re = input('Choose maximum value of Re: ');
Step = input('Choose length of steps of Re: ');
Re = Min_Re:Step:Max_Re;

This creates a vector for Re (based on user given interval and user given step value) and I need this so as to solve the equation for every value of Re and then when having done this create a plot(Re, f) showing graphically how f changes with Re.

So far I've used fzero without any success.

This is how I wrote my function:

function y = y(f)
Min_Re = input('Choose minimum value of Re: ');
Max_Re = input('Choose maximum value of Re: ');
Step = input('Choose length of step: ');
Re = Min_Re:Step:Max_Re;

y = 4*log10(Re.*sqrt(f)) - 0.4 - 1/sqrt(f);

Then the program/script actually trying to solve the equation using fzero as:

f0 = input('Choose starting guess: ');

f = fzero(@(f) y(f), f0);

This is the error message I get:

Operands to the || and && operators must be convertible to logical scalar values.

Error in fzero (line 308)
    elseif ~isfinite(fx) || ~isreal(fx)

Error in untitled2 (line 3)
f = fzero(@(f) y(f), f0);

I really need to solve this problem but most importantly I do not want any new and complex built-in functions that I must learn/understand and then also be able to explain.

This is by the way the message I get from MATLAB when I try using only one value for Re (Re = 10000):

Choose starting guess: 10
Exiting fzero: aborting search for an interval containing a sign change
    because complex function value encountered during search.
(Function value at -2.8 is 16.4943+3.32637i.)
Check function or try again with a different starting value.

Appreciate all the help I can get.

Thank you for your time and thank you in advance.

Cheers!

[Matt J]

Please put your code and Error messages in a separate font from your text using this toolbar button

Answer 1

[Roger Stafford]

It may be of some interest to you that solutions to your equation can be expressed explicitly in terms of the 'lambertw' function. Starting with your original equation it can be transformed from

 1/sqrt(f)-4*log10(Re*sqrt(f))-.4 = 0

to the form

 1/sqrt(f)-k1-k2*log(sqrt(f)) = 0

using k1 = 4/log(10)*log(Re)+.4 and k2 = 4/log(10). Then substituting log(sqrt(f)) = x which is equivalent to sqrt(f) = exp(x) and to f = exp(2*x), it takes on the form

 1/exp(x)-k1-k2*x = 0

Your Symbolic Toolbox is then capable of solving this equation for x explicitly:

 x = lambertw(exp(k1/k2)/k2)-k1/k2

I believe you can find the 'lambertw' function somewhere in either the Symbolic Toolbox or in Mupad. It has the characteristic that for real solutions there is one range of arguments which give single values, another range for which there are two branches (two values) for each argument, and a third range for which there are no solutions. Hopefully for all your values of the Re in k1 you will have solutions. From such solutions you can work backwards to the corresponding values of f using the above relationship,

 f = exp(2*x).

[Brod]

Thank you for you help but I feel that this is not the way that we are supposed to solve this problem. This is part of course in Scientific Computing and this part of the course is especially for solving non-linear equations numerically.

But I will take a closer look at your suggesting.

Also, I do not quite follow how you transformed by original equation. Especially the part with using k1 = 4/log(10) (this is isn't log(10) this is the 10-logarithm that is in the equation.

[Roger Stafford]

That follows from the fact that

 log10(x) = log(x)/log(10)

with the understanding that in matlab 'log' means natural logarithm.

The trouble you had with 'fzero' was due to trying call it to do many things with just one call. 'fzero' won't work like that. You have to use a loop to call it separately for each number it solves, that is in your case at each 'step'.

[Brod]

So how would you propose creating a loop to do this. I figured letting it do like this:

a = 1;

while a < Max_Re
a = a + 1
end

But should I put this in the script of the function? Or in the equation solving function?

If it is in the script of the function then I would change y to:

y = 4*log10(Re(a).*sqrt(f)) - 0.4 - 1/sqrt(f);

And also, how would I create the very much needed vector for f?