|
yup Newton Raphsom method would do it.
Also It would be helpful to find out what are the ranges of possible values of phi, say if phi belongs in -1 to 1 you can apply sine series formula to phi.
|
|
|
|
|
The equation is modeling a physical phenomenon (water streaming in a pipe), I don't know all the details, however I'm confident it will be well behaved, have a simple first guess (say 10 degrees), and converge easily.
|
|
|
|
|
What do the variables represent? Is there a way to determine reasonable variable limits? For example, if you know that Φ must be between +/- 2 Π that limits the solutions and table look up or a curve fit to the solution might be good enough. How much precision do you need?
CQ de W5ALT
Walt Fair, Jr., P. E.
Comport Computing
Specializing in Technical Engineering Software
|
|
|
|
|
I was thinking perhaps a cubic spline solution...
The variables are part of the equation for fluid flow in a pipe, and Φ represents the angle from horizontal a line drawn from the center of the pipe to the intersection of the line marking the top of the fluid within.
V = (1.485/n)* R2/3 * S1/2, ft/sec.
S is the slope of the pipe, and n is a very small number that varies with the smoothness of the pipe - a common value is .013. R is called the 'hydraulic radius' and is given by the ratio of the area of flow to the wetted perimeter in the pipe. Calculating that area value is what hurts my head.
It's fairly straightforward to calculate the flow, given everything else, but I'm trying to predict a level in the pipe from a given flow and slope. That's not what the equation is designed for, though there's no reason not to use it this way. If I can't find a reasonable way to do it, I may have to resort to Bernoulli's equation, with all those nasty energy calculations and differential equations. We hates those...
"A Journey of a Thousand Rest Stops Begins with a Single Movement"
|
|
|
|
|
Roger Wright wrote: I was thinking perhaps a cubic spline solution...
A cubic spline would be an approximation, i.e. a curve fitting a couple of points. Calculating it would be more complex than solving the equation itself, and the results would not be all that accurate.
|
|
|
|
|
Sure, Manning's equation. It's been awhile since I worked with single phase open channel flow. S should be the slope of the hydraulic head (potential gradient) if I recall. That's not the same as the pipe slope, unless the head loss is negligible or the pipe is full. If it's full, then Darcy-Weisbach is a better way to estimate and if you're in fully turbulent flow, the friction factors can be estimated quite easily, but it's still probably not an analytic solution.
I'm not sure I understand what Φ is in relation to the hydraulic gradient and/or pipe slope, though. Where does the Φ come in?
CQ de W5ALT
Walt Fair, Jr., P. E.
Comport Computing
Specializing in Technical Engineering Software
|
|
|
|
|
I doubt that Manning's is the right equation to use, but having searched out a number of public building codes I've found that all of them require sewer flows to be calculated using it. Even crazier, most of them specify that n = .013 be used, regardless of the pipe material. Darcy-Weisbach is probably better, but the only application I've found using it requires reading values from a nomograph - ridiculous! IIRC, I believe it requires knowledge of the viscosity and flow characteristics, as well, little of which is known at design time. I can estimate %solids from the measurements taken from the influent at the sewer processing plant, but even then I probably have areas with wildly different fluid properties.
The angle Φ is used to calculate the area of fluid in a pipe that isn't full. Picture a circle with a line drawn through it representing the fluid level. A line from the center of the circle to the point where the fluid level intersects the circle makes the angle Φ with the horizontal.
I've got some more figuring to do, I guess. Or perhaps I should just keep it simple.
"A Journey of a Thousand Rest Stops Begins with a Single Movement"
|
|
|
|
|
OK, I understand, Φ is based on Θ in this figure where they calculate hydraulic radius: http://www.ajdesigner.com/phphydraulicradius/hydraulic_radius_equation_pipe.php[^]. If I follow you, Φ = (Π - Θ)/2.
Since it's building codes that are specifying the method to use, then you have little choice but to do that, unless you have a licensed engineer come up with a good reason to do something better. Manning's equation with the empirical friction factor equivalent is probably conservative for the purpose, and that's OK.
I'd guess that the iterative solution is your best bet.
CQ de W5ALT
Walt Fair, Jr., P. E.
Comport Computing
Specializing in Technical Engineering Software
|
|
|
|
|
Walt Fair, Jr. wrote: I'd guess that the iterative solution is your best bet.
we are all converging now. 
|
|
|
|
|
I made this, but it's not an analytical answer..
static double Solve(double y)
{
double x = y;
double error = double.MaxValue;
double prevError = double.PositiveInfinity;
while (error != 0 && Math.Abs(prevError) > Math.Abs(error))
{
x -= (x - Math.Sin(x) - y) / (1 - Math.Cos(x));
prevError = error;
error = y - (x - Math.Sin(x));
}
return x;
}
It doesn't work for zero, and gives infinity for y < 0.00000001054, but otherwise it seems alright.
The input is supposed to be your value for 4R/D and the output will approximate Φ.
I can't guarantee that it will always terminate, but it did for every value I tested (which included negative numbers, and insanely big numbers, and very small numbers)
|
|
|
|
|
do I detect a Newton-Raphson smell around here?
|
|
|
|
|
Probably, because that's what it is
|
|
|
|
|
Thanks! That should be a good start. I may have a shortcut to make it converge faster, though. I'll have to think about it.
The idea is to model (one situation) a sewer line that is existing, to which a new customer - a car wash, for instance - wants to connect. I'd like to then force a change in gpm inflow, then see what happens to everything else. Since the gpm is a function of the area and velocity of flow, that's one equation. Then the other can be used to provide a bound for the next guess at Φ until the two equations converge to a single value, +/- some tolerance value.
My head hurts; I think I'll bang out something easy - a database program to inventory the entire electrical system maybe - while I'm cogitating. 
"A Journey of a Thousand Rest Stops Begins with a Single Movement"
|
|
|
|
|
That went right over my head..
Maybe I'll think about it some more but I have other things to do as well.
Good luck 
|
|
|
|
|
Before you start optimizing things, may I suggest you observe the current behavior. In my experience, equations that model macroscopic physical phenomena are well-behaved and tend to converge quite well. So just log each iteration, and you will probably notice a few iterations will bring your error down to less than 0.1%, which is probably much better than the accuracy of the formula anyway.
|
|
|
|
|
x - sin(x) = y
x = pi/2*(y/(pi/2-1)^(1/3)
accurate to 3.8% for 0 <= y <= pi/2-1
follow this by one NR refinement (for y > 0)
x1 = x0 + (y - x0 + sin(x0))/(1 - cos(x0))
reduces relative error to less than 1.5 x 10^-3 over this range
two refinements gets relative error down to 2 x 10^-6
three refinements gets relative error down to 5 x 10^-12
Peter
"Until the invention of the computer, the machine gun was the device that enabled humans to make the most mistakes in the smallest amount of time."
|
|
|
|
|
That's beautiful, Peter. Thank you! God bless math wizards! 
[EDIT]
A question:
The term (2-1) appears often. That makes no sense, so I'm guessing that it's an HTML encoding error. Can you clarify, please?
"A Journey of a Thousand Rest Stops Begins with a Single Movement"
modified on Thursday, June 24, 2010 2:56 AM
|
|
|
|
|
Actually a slightly more inspired approximation is
x - sin(x) = y
let q = pi/2
calculate b = (q^3/(q-1)-6)/(q-1) approx equal to 1.384269281
then
x = ((6 + b*y) * y)^(1/3)
is accurate to 0.47% for 0 <= y <= pi/2-1
follow this by one NR refinement (for y > 0)
x1 = x0 + (y - x0 + sin(x0))/(1 - cos(x0))
reduces relative error to less than 2 x 10^-5 over this range
two refinements gets relative error down to 4 x 10^-10
re the 2-1, just evaluate is as your compiler would, so pi/2-1 = (pi/2)-1
Peter
"Until the invention of the computer, the machine gun was the device that enabled humans to make the most mistakes in the smallest amount of time."
|
|
|
|
|
How are you making these approximations?
|
|
|
|
|
harold aptroot wrote: How are you making these approximations?
for x - sin(x) = y, for small y
x = (6y)^(1/3)
the first approximation is simply taking this form and making the value correct at the two ends of the interval, y in (0, ym) where ym = pi/2-1. x=0 for y=0, so we only have to make x=pi/2 at y=ym, so
x = pi/2 * (y/ym)^(1/3)
This gives about 4% error.
Noting that the scale factor is wrong at y=0, the second approximation is of the form
x = ((6+b*y)*y)^(1/3)
which is like x = (6y)^(1/3) for small y, and we choose b so that x=pi/2 for y=ym. This gives less than 0.5% error.
Add this to a simple excel spreadsheet to do a few NR iterations and calculate the errors and you are done.
I'm sure that you could improve things given enough time.
Peter
"Until the invention of the computer, the machine gun was the device that enabled humans to make the most mistakes in the smallest amount of time."
|
|
|
|
|
cp9876 wrote: x = ((6+b*y)*y)^(1/3)
which is like x = (6y)^(1/3) for small y, and we choose b so that x=pi/2 for y=ym. This gives less than 0.5% error.
I'm sure that you could improve things given enough time.
Using the same b and ym as before, put
x = ((6+b*ym*(y/ym)^a)*y)^(1/3)
Run Excel's solver to optimise 'a' gives a = 0.731889437 for a relative error < 0.033%
One NR iteration reduces this to 1 x 10^-7
OK, I'm going to stop now....
Peter
"Until the invention of the computer, the machine gun was the device that enabled humans to make the most mistakes in the smallest amount of time."
|
|
|
|
|
How did you know what "form" the formula should have?
|
|
|
|
|
harold aptroot wrote: How did you know what "form" the formula should have?
if you are asking where the
x = (6y)^(1/3)
comes from, it is just a first order expansion around y=0
Just to summarise for Roger, the best approxmiation so far:
let q = pi/2
calculate b1 = (q^3/(q-1)-6)
ym = q-1
x = ((6+b1*(y/ym)^a)*y)^(1/3)
with a = 0.731889437
gives a relative error < 0.033%
One NR iteration reduces this to 1 x 10^-7
Peter
"Until the invention of the computer, the machine gun was the device that enabled humans to make the most mistakes in the smallest amount of time."
|
|
|
|
|
Excellent! I have no clear idea about how it is derived, but that won't stop me from using what works.
When you refer to a first order expansion, what is it you're expanding? Taylor series?
"A Journey of a Thousand Rest Stops Begins with a Single Movement"
|
|
|
|
|
Roger Wright wrote: When you refer to a first order expansion, what is it you're expanding? Taylor series?
Just do a Taylor expansion
x - sin(x) = x - (x - x^3/3! + ..) = y
as x->0 only the lowest powers are significant, so
x^3 / 6 = y
or x = (6y)^(1/3)
The rest of the approximation stuff is guesswork, and using something like Excel to see how well you are going. It's a bit of a game that's fun to play occasionally.
Using Excel and the solve function you can get a nice tidy approximation that's probably better for implementation:
x = (ay)^(1/3) + by + cy^2
a = 5.998854
b = 0.10235
c = 0.01579
the relative error is less than 2.7 x 10^-5 over the range 0 <= y <= pi/2 - 1
If you need better than this, one NR iteration gets you to 10^-10 territory and one more should pretty well get you full double precision.
Peter
"Until the invention of the computer, the machine gun was the device that enabled humans to make the most mistakes in the smallest amount of time."
|
|
|
|
General 
News 
Suggestion 
Question 
Bug 
Answer 
Joke 
Praise 
Rant 
Admin 
Submit an article or tip
