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_med=Q_diario/(24*3600) %dividimos el caudal a lo largo del dia y lo pasamos a m3/s
Q_punta=coef_punta*Q_med
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)
%%
%dimensiones de tubería a aliviar
%tomo tubería rectangular para este ejemplo
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
n=0.015
%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)
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 offIm 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_med=Q_diario/(24*3600) %dividimos el caudal a lo largo del dia y lo pasamos a m3/s
Q_punta=coef_punta*Q_med
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)
%%
%dimensiones de tubería a aliviar
%tomo tubería rectangular para este ejemplo
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
n=0.015
%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)
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 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_med=Q_diario/(24*3600) %dividimos el caudal a lo largo del dia y lo pasamos a m3/s
Q_punta=coef_punta*Q_med
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)
%%
%dimensiones de tubería a aliviar
%tomo tubería rectangular para este ejemplo
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
n=0.015
%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)
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 fzero, convergence, exitflag=1 MATLAB Answers — New Questions
