"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.866569319.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 yvalue 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 CivE!
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
%%%%%
