Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
[Solve] Calculating uniform depth with strickler

Subject: [Solve] Calculating uniform depth with strickler

From: Maarten

Date: 10 Dec, 2010 17:19:07

Message: 1 of 5

Hello all,

For a civil engineering student research to a new bypass river for high water somewhere in the Netherlands, we have to calculate uniform depth and length of the backwater curve.
Normally that would not be a problem, but we want to do it fully automated, because we have to make a lot of calculations by hand.

In general: these are the formulae given for the calculation:
A = (b+m*y)*y (wet surface)
R = A/(b+2*y*sqrt(1+m^2)) (hydraulic radius)
Q = k*A*R^(2/3)*s^(1/2) (strickler formula)

where
y = water level (which we have to calculate)
b = width of the river bed [m]
m = slope of the dikes
A = wet surface [m2]
R = hydraulic radius [m]
k = strickler coëfficiënt [m^(1/3)/s]
s = slope of river
Q = flow rate of the river [m3/s]

With given:
s=0.0001
k=25
b=360
m=2

Calculating y by hand it gives us an uniform depth of 18.55 meters, this is about realistic.

Calculting y with Maple 13 using solve gives us an irreal number which is -186.8665693-19.10925842*i.
However, when using fsolve instead of solve in Maple 13 this does give us 18.55 again.
A simple plot (Q as a function of y with given data) shows this is correct

However, now I'm trying to do that in Matlab as well. It's not going too well.

So with many effort I created 4 matlab files, because I need functions for A, R and the strickler formula.
Because strickler uses R and A en R uses A I had to triple define A and double define R. Is there an easier way (by for example crossrefering to different files, I need all the functions in the end)?

However.

A.m

function z=A(y)
    global m;
    global b;
    z=(b+m*y)*y;
end

R.m

function z=R(y)
global b;
global m;
function z=A(y)
      z=(b+m*y)*y;
end
z=A(y)/(b+2*y*sqrt(1+m^2));
end

strickler.m

function z=strickler(y)
global b;
global m;
global s;
global k;
    function z=A(y)
    z=(b+m*y)*y;
    end
    function z=R(y)
    z=A(y)/(b+2*y*sqrt(1+m^2));
    end
z=k*A(y)*R(y)^(2/3)*s^(1/2);
end

And the last file used as the model:

werkendam.m

global m;
global b;
global s;
global k;

m=2;
b=360;
s=0.0001;
k=25;

Q=12000;

syms y;
solve('stricler(y)-12000',y)

Executing this latest file to solve the y-value it gives a rather complex error:

??? Error using ==> mupadmex
Error in MuPAD command: cannot differentiate equation [numeric::fsolve]

Error in ==> sym.sym>sym.mupadmexnout at 2018
        out = mupadmex(fcn,args{:});

Error in ==> solve at 76
[symvars,R] = mupadmexnout('symobj::solvefull',eqns,vars);

Error in ==> werkendam at 16
solve('stricler(y)-12000',y)

When replacing in solve('stricler(y)-12000',y) the 12000 by Q, resulting in:
solve('stricler(y)-Q',y), gives me as an output:

Warning: Explicit solution could not be found.
> In solve at 81
  In werkendam_test at 16
 
ans =
 
[ empty sym ]

I just cannot figure out what's wrong!
Could somebody please help me?
Many thanks!

Subject: [Solve] Calculating uniform depth with strickler

From: Sean de

Date: 10 Dec, 2010 17:56:05

Message: 2 of 5

"Maarten " <m.a.schoemaker@student.tudelft.nl> wrote in message <idtneb$fso$1@fred.mathworks.com>...
> Hello all,
>
> For a civil engineering student research to a new bypass river for high water somewhere in the Netherlands, we have to calculate uniform depth and length of the backwater curve.
> Normally that would not be a problem, but we want to do it fully automated, because we have to make a lot of calculations by hand.
>
> In general: these are the formulae given for the calculation:
> A = (b+m*y)*y (wet surface)
> R = A/(b+2*y*sqrt(1+m^2)) (hydraulic radius)
> Q = k*A*R^(2/3)*s^(1/2) (strickler formula)
>
> where
> y = water level (which we have to calculate)
> b = width of the river bed [m]
> m = slope of the dikes
> A = wet surface [m2]
> R = hydraulic radius [m]
> k = strickler coëfficiënt [m^(1/3)/s]
> s = slope of river
> Q = flow rate of the river [m3/s]
>
> With given:
> s=0.0001
> k=25
> b=360
> m=2
>
> Calculating y by hand it gives us an uniform depth of 18.55 meters, this is about realistic.
>
> Calculting y with Maple 13 using solve gives us an irreal number which is -186.8665693-19.10925842*i.
> However, when using fsolve instead of solve in Maple 13 this does give us 18.55 again.
> A simple plot (Q as a function of y with given data) shows this is correct
>
> However, now I'm trying to do that in Matlab as well. It's not going too well.
>
> So with many effort I created 4 matlab files, because I need functions for A, R and the strickler formula.
> Because strickler uses R and A en R uses A I had to triple define A and double define R. Is there an easier way (by for example crossrefering to different files, I need all the functions in the end)?
>
> However.
>
> A.m
>
> function z=A(y)
> global m;
> global b;
> z=(b+m*y)*y;
> end
>
> R.m
>
> function z=R(y)
> global b;
> global m;
> function z=A(y)
> z=(b+m*y)*y;
> end
> z=A(y)/(b+2*y*sqrt(1+m^2));
> end
>
> strickler.m
>
> function z=strickler(y)
> global b;
> global m;
> global s;
> global k;
> function z=A(y)
> z=(b+m*y)*y;
> end
> function z=R(y)
> z=A(y)/(b+2*y*sqrt(1+m^2));
> end
> z=k*A(y)*R(y)^(2/3)*s^(1/2);
> end
>
> And the last file used as the model:
>
> werkendam.m
>
> global m;
> global b;
> global s;
> global k;
>
> m=2;
> b=360;
> s=0.0001;
> k=25;
>
> Q=12000;
>
> syms y;
> solve('stricler(y)-12000',y)
>
> Executing this latest file to solve the y-value it gives a rather complex error:
>
> ??? Error using ==> mupadmex
> Error in MuPAD command: cannot differentiate equation [numeric::fsolve]
>
> Error in ==> sym.sym>sym.mupadmexnout at 2018
> out = mupadmex(fcn,args{:});
>
> Error in ==> solve at 76
> [symvars,R] = mupadmexnout('symobj::solvefull',eqns,vars);
>
> Error in ==> werkendam at 16
> solve('stricler(y)-12000',y)
>
> When replacing in solve('stricler(y)-12000',y) the 12000 by Q, resulting in:
> solve('stricler(y)-Q',y), gives me as an output:
>
> Warning: Explicit solution could not be found.
> > In solve at 81
> In werkendam_test at 16
>
> ans =
>
> [ empty sym ]
>
> I just cannot figure out what's wrong!
> Could somebody please help me?
> Many thanks!


Fellow Civ-E!

The reason it could not solve was because if y = 0, the whole expression = 0. Thus instead of running the engine in symbolic notation, run it in vpa with a high initial guess. You could also do this whole thing strictly numerically without symbolic variables, though this works fine.

%%%%%
syms b m y s k w

A = (b+m*y)*y;
R = A/(b+2*y*sqrt(1+m^2));
Q = k*A*R^(2/3)*s^(1/2);

Qs = vpa(subs(Q,{m b s k},{2 360 .0001 25})); %Expression
fzero(@(y)subs(Qs,'y',y)-12000,20) %Solve it with an initial guess of 20

%SCd
%%%%%

Subject: [Solve] Calculating uniform depth with strickler

From: Maarten

Date: 10 Dec, 2010 19:36:05

Message: 3 of 5

"Sean de " <sean.dewolski@nospamplease.umit.maine.edu> wrote in message <idtpjl$nqo$1@fred.mathworks.com>...
>
> Fellow Civ-E!
>
> The reason it could not solve was because if y = 0, the whole expression = 0. Thus instead of running the engine in symbolic notation, run it in vpa with a high initial guess. You could also do this whole thing strictly numerically without symbolic variables, though this works fine.
>
> %%%%%
> syms b m y s k w
>
> A = (b+m*y)*y;
> R = A/(b+2*y*sqrt(1+m^2));
> Q = k*A*R^(2/3)*s^(1/2);
>
> Qs = vpa(subs(Q,{m b s k},{2 360 .0001 25})); %Expression
> fzero(@(y)subs(Qs,'y',y)-12000,20) %Solve it with an initial guess of 20
>
> %SCd
> %%%%%

ans =

   18.5519

Many, many thanks!

Part of this research (3rd year student) is also learning these things.
(BTW: greetings from TU Delft)

Subject: [Solve] Calculating uniform depth with strickler

From: Maarten

Date: 10 Dec, 2010 19:44:20

Message: 4 of 5

"Sean de " <sean.dewolski@nospamplease.umit.maine.edu> wrote in message <idtpjl$nqo$1@fred.mathworks.com>...
>
> Fellow Civ-E!
>
> The reason it could not solve was because if y = 0, the whole expression = 0. Thus instead of running the engine in symbolic notation, run it in vpa with a high initial guess. You could also do this whole thing strictly numerically without symbolic variables, though this works fine.
>
> %%%%%
> syms b m y s k w
>
> A = (b+m*y)*y;
> R = A/(b+2*y*sqrt(1+m^2));
> Q = k*A*R^(2/3)*s^(1/2);
>
> Qs = vpa(subs(Q,{m b s k},{2 360 .0001 25})); %Expression
> fzero(@(y)subs(Qs,'y',y)-12000,20) %Solve it with an initial guess of 20
>
> %SCd
> %%%%%

I sort of adapted it as well in order to be able to change the variables easier:

syms y

b=360;
s=0.0001;
k=25;
m=2;

A = (b+m*y)*y;
R = A/(b+2*y*sqrt(1+m^2));
Q = k*A*R^(2/3)*s^(1/2);

fzero(@(y)subs(Q,'y',y)-12000,20) %Solve it with an initial guess of 20

Still works great!

Subject: [Solve] Calculating uniform depth with strickler

From: Steven_Lord

Date: 10 Dec, 2010 23:01:35

Message: 5 of 5



"Maarten " <m.a.schoemaker@student.tudelft.nl> wrote in message
news:idtvuk$kh7$1@fred.mathworks.com...
> "Sean de " <sean.dewolski@nospamplease.umit.maine.edu> wrote in message
> <idtpjl$nqo$1@fred.mathworks.com>...

*snip*

> I sort of adapted it as well in order to be able to change the variables
> easier:
>
> syms y
>
> b=360;
> s=0.0001;
> k=25;
> m=2;
>
> A = (b+m*y)*y;
> R = A/(b+2*y*sqrt(1+m^2));
> Q = k*A*R^(2/3)*s^(1/2);
>
> fzero(@(y)subs(Q,'y',y)-12000,20) %Solve it with an initial guess of 20
>
> Still works great!

You might be interested in the MATLABFUNCTION function to directly create an
anonymous function handle for use by FZERO from your symbolic expression.

--
Steve Lord
slord@mathworks.com
comp.soft-sys.matlab (CSSM) FAQ: http://matlab.wikia.com/wiki/FAQ
To contact Technical Support use the Contact Us link on
http://www.mathworks.com

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us