Fzero not converging, giving previous answer back.

Im trying to solve the water level across a fixed lip spillway by using fzero, but I ran into the problem that the function gives an initial answer and doesn't change between iterations in my for loop. By manually advancing throuhg the loops i can see that in the first one fzero already begins drifting and stops converging on an answer, giving me my last iterations (the first ever answer it got) back to me.

close all
clear all 
clc
%% Que longitud de aliviadero necesito para aliviar hasta un calado determinado con la entrada de un caudal dado?

%%
%determino caudal de entrada
Poblacion=10000;
Dotacion=200; %l/hab*dia
coef_aguas_negras=0.8;%cuanto agua dotada pasa a aguias negras
coef_punta=2; %coeficiente de punta de aguas residuales pico comparada con caudal medio, 2 en poblaciones pequeñas (desfavorable)
Q_diario=Poblacion*Dotacion*coef_aguas_negras*1/1000%pasamos a m3/dia

Q_diario = 1600

Q_med=Q_diario/(24*3600) %dividimos el caudal a lo largo del dia y lo pasamos a m3/s

Q_med = 0.0185

Q_punta=coef_punta*Q_med

Q_punta = 0.0370

Q_gris=10*Q_med%caudal a partir del cual se determina agua gris y está suficientemente diluido como para verter a cauce sin tratar (aliviadero vierte a partir de este caudal)

Q_gris = 0.1852

%%
%dimensiones de tubería a aliviar
%tomo tubería rectangular para este ejemplo
b=1

b = 1

%para calcular c se necesita calcular la velocidad de flujo por el conducto
%(V=Q/A)
i=0.001 %pendiente del fondo del 1 por mil

i = 1.0000e-03

n=0.015

n = 0.0150

%uso manning para conseguir c
ec1=@(y) 1/n*((b*y)/(b+2*y))^(2/3)*i^(1/2)*b*y-Q_gris;
y=fzero(ec1,1);
%redondeo hacia arriba (los 10 cm hacia arriba)
c=ceil(y*10)/10;
%%
%creo el espacio del vector de puntos diferenciales a lo largo del
%aliviadero
L=5;
dx=0.1;
x=[0:dx:L];
%%
%INPUT
%meto el nivel al principio del aliviadero (se saca igual con manning
%arriba si tenemos un caudal de entrada
Q=1;
ec2=@(y) 1/n*((b*y)/(b+2*y))^(2/3)*i^(1/2)*[b*y]-Q;
%%
h1=zeros(1,length(x)-1);
h2=zeros(1,length(x)-1);
v1=zeros(1,length(x)-1);
v2=zeros(1,length(x)-1);

y1=fzero(ec2,1);
h1(1)=y1;
v_e=Q/(b*h1(1));
v1(1)=v_e;
q1=Q;
Q1(1)=q1;
Cv=1.86; %para la formula de aliviadero de pared delgada
g=9.81;
ec3=@(y2) 10^16*([y1]+[((v_e)^2)/(2*g)]-[y2]-[(((q1-Cv*dx*(((y1+y2)/2)-c)^(3/2))/(b*y2))^2)/(2*g)]-[sqrt(((((v_e+((q1-Cv*dx*(((y1+y2)/2)-c)^(3/2))/(b*y2))))/2)*n*((b+2*y2)^(2/3)))/((b*y2)^(2/3)))*dx]);

a=fzero(ec3,10);
y2=a;
h2(1)=y2;
q2=(q1-Cv*dx*(((y1+y2)/2)-c)^(3/2));
Q2(1)=q2;
v2(1)=q2/(b*h2(1));
%%
for k=2:length(x)
    y1=h2(k-1);
    v_e=v2(k-1);
    q1=Q2(k-1);
    a=fzero(ec3,10);
    y2=a;
    h2(k)=y2;
    q2=(q1-Cv*dx*(((y1+y2)/2)-c)^(3/2));
    Q2(k)=q2;
    v2(k)=q2/(b*h2(k));
    h1(k)=y1;
    v1(k)=v_e;
end
%%
%intento diagnosticar el error de convergencia de fzero

figure ;
fplot(ec3,[-2,2]);

Warning: Function behaves unexpectedly on array inputs. To improve performance, properly vectorize your function to return an output with the same size and shape as the input arguments.

hold on
plot(y2,0,'r*');
yline(0);
hold off

k=2;
y1=h2(k-1);
    v_e=v2(k-1);
    q1=Q2(k-1);
    [a,fval,exitflag,output]=fzero(ec3,10)

a = 0.9815

fval = 0.5204

exitflag = 1

output = struct with fields:
    intervaliterations: 11
            iterations: 6
             funcCount: 28
             algorithm: 'bisection, interpolation'
               message: 'Zero found in the interval [0.949033, 16.4]'

    drift=[y1]+[((v_e)^2)/(2*g)]-[y2]-[(((q1-Cv*dx*(((y1+y2)/2)-c)^(3/2))/(b*y2))^2)/(2*g)]-[sqrt(((((v_e+((q1-Cv*dx*(((y1+y2)/2)-c)^(3/2))/(b*y2))))/2)*n*((b+2*y2)^(2/3)))/((b*y2)^(2/3)))*dx];

    figure ;
fplot(ec3,[-2,2]);

Warning: Function behaves unexpectedly on array inputs. To improve performance, properly vectorize your function to return an output with the same size and shape as the input arguments.

hold on
plot(y2,0,'r*');
yline(0);
hold off

[Ivo]

The problem is with fzero solving ec3, i checked the equation and it's correctly writen, but by checking the value of 'drift' you can see that the answers it gives do not zero the equation. With the exitflag instruction it also tells me it doesn't converge on an answer.

The original ec3 did not have the 10^16, but I added that after seeing that maybe the problem was with the tolerance (the fval for the answers where below the 2.22e^-16 of tolerance, so i thought it was that), but even scaling the values up doesnt seem to change the answer (fval is now in the single digits range and still doesnt converge).

[Walter Roberson]

Fixing the fplot() warnings

close all
clear all 
clc
%% Que longitud de aliviadero necesito para aliviar hasta un calado determinado con la entrada de un caudal dado?

%%
%determino caudal de entrada
Poblacion=10000;
Dotacion=200; %l/hab*dia
coef_aguas_negras=0.8;%cuanto agua dotada pasa a aguias negras
coef_punta=2; %coeficiente de punta de aguas residuales pico comparada con caudal medio, 2 en poblaciones pequeñas (desfavorable)
Q_diario=Poblacion*Dotacion*coef_aguas_negras*1/1000%pasamos a m3/dia

Q_diario = 1600

Q_med=Q_diario/(24*3600) %dividimos el caudal a lo largo del dia y lo pasamos a m3/s

Q_med = 0.0185

Q_punta=coef_punta*Q_med

Q_punta = 0.0370

Q_gris=10*Q_med%caudal a partir del cual se determina agua gris y está suficientemente diluido como para verter a cauce sin tratar (aliviadero vierte a partir de este caudal)

Q_gris = 0.1852

%%
%dimensiones de tubería a aliviar
%tomo tubería rectangular para este ejemplo
b=1

b = 1

%para calcular c se necesita calcular la velocidad de flujo por el conducto
%(V=Q/A)
i=0.001 %pendiente del fondo del 1 por mil

i = 1.0000e-03

n=0.015

n = 0.0150

%uso manning para conseguir c
ec1=@(y) 1/n*((b*y)./(b+2*y)).^(2/3)*i^(1/2)*b.*y-Q_gris;
y=fzero(ec1,1);
%redondeo hacia arriba (los 10 cm hacia arriba)
c=ceil(y*10)/10;
%%
%creo el espacio del vector de puntos diferenciales a lo largo del
%aliviadero
L=5;
dx=0.1;
x=[0:dx:L];
%%
%INPUT
%meto el nivel al principio del aliviadero (se saca igual con manning
%arriba si tenemos un caudal de entrada
Q=1;
ec2=@(y) 1/n*((b*y)./(b+2*y)).^(2/3)*i^(1/2).*[b.*y]-Q;
%%
h1=zeros(1,length(x)-1);
h2=zeros(1,length(x)-1);
v1=zeros(1,length(x)-1);
v2=zeros(1,length(x)-1);

y1=fzero(ec2,1);
h1(1)=y1;
v_e=Q/(b*h1(1));
v1(1)=v_e;
q1=Q;
Q1(1)=q1;
Cv=1.86; %para la formula de aliviadero de pared delgada
g=9.81;
ec3=@(y2) 10^16*([y1]+[((v_e)^2)/(2*g)]-[y2]-[(((q1-Cv*dx.*(((y1+y2)/2)-c).^(3/2))./(b*y2)).^2)/(2*g)]-[sqrt(((((v_e+((q1-Cv*dx*(((y1+y2)/2)-c).^(3/2))./(b*y2))))/2)*n.*((b+2*y2).^(2/3)))./((b*y2).^(2/3)))*dx]);

a=fzero(ec3,10);
y2=a;
h2(1)=y2;
q2=(q1-Cv*dx*(((y1+y2)/2)-c)^(3/2));
Q2(1)=q2;
v2(1)=q2/(b*h2(1));
%%
for k=2:length(x)
    y1=h2(k-1);
    v_e=v2(k-1);
    q1=Q2(k-1);
    a=fzero(ec3,10);
    y2=a;
    h2(k)=y2;
    q2=(q1-Cv*dx*(((y1+y2)/2)-c)^(3/2));
    Q2(k)=q2;
    v2(k)=q2/(b*h2(k));
    h1(k)=y1;
    v1(k)=v_e;
end
%%
%intento diagnosticar el error de convergencia de fzero

figure ;
fplot(ec3,[-2,2]);
hold on
plot(y2,0,'r*');
yline(0);
hold off

●

k=2;
y1=h2(k-1);
    v_e=v2(k-1);
    q1=Q2(k-1);
    [a,fval,exitflag,output]=fzero(ec3,10)

a = 0.9815

fval = 0.5204

exitflag = 1

output = struct with fields:
    intervaliterations: 11
            iterations: 6
             funcCount: 28
             algorithm: 'bisection, interpolation'
               message: 'Zero found in the interval [0.949033, 16.4]'

    drift=[y1]+[((v_e)^2)/(2*g)]-[y2]-[(((q1-Cv*dx*(((y1+y2)/2)-c)^(3/2))/(b*y2))^2)/(2*g)]-[sqrt(((((v_e+((q1-Cv*dx*(((y1+y2)/2)-c)^(3/2))/(b*y2))))/2)*n*((b+2*y2)^(2/3)))/((b*y2)^(2/3)))*dx];

    figure ;
fplot(ec3,[-2,2]);
hold on
plot(y2,0,'r*');
yline(0);
hold off

Answer 1

[Torsten]

The constants used in your function "ec3" are not automatically changed if you recompute them in your code. So each time before you use "ec3" in the call to "fsolve", you must redefine the function "ec3". This will ensure that the constants in the function are set to their actual values.

[Ivo]

Oh, that would explain why it only gives one result, thanks. I'll try that and see if it works.