The derivative needed for h1
The "trick" to an Ag Engineering problem was in solving for h1.
Given:
dh=3
bw=12
cr=4
cs=0.01
ss1=6
ss2=6
ss3=6
CF=1.5
ac=((bw^2*cs+2*bw*(dh-h1)*(1+cs*ss2)+(dh-h1)^2*(1+cs*ss2)*(ss2+ss3))/(2-2*cs*ss3))
ad=((cr^2*cs+h1^2*(ss1+ss2)*(1+cs*ss2)+2*cr*(h1+cs*h1*ss2))/(2-2*cs*ss1))
CF=ac/ad
Find h1: h1 is about 1.5.
There are several ways to find the true value of h1. I showed one in the solution to "An Ag Engineering problem".
But a good way would need the derivative of how CF changes with h1. I finally solved that. And with Newton's
method (also known as the Newton–Raphson method, the solution is quick and accurate.
Can anyone show the derivative or show a better method?
Given:
dh=3
bw=12
cr=4
cs=0.01
ss1=6
ss2=6
ss3=6
CF=1.5
ac=((bw^2*cs+2*bw*(dh-h1)*(1+cs*ss2)+(dh-h1)^2*(1+cs*ss2)*(ss2+ss3))/(2-2*cs*ss3))
ad=((cr^2*cs+h1^2*(ss1+ss2)*(1+cs*ss2)+2*cr*(h1+cs*h1*ss2))/(2-2*cs*ss1))
CF=ac/ad
Find h1: h1 is about 1.5.
There are several ways to find the true value of h1. I showed one in the solution to "An Ag Engineering problem".
But a good way would need the derivative of how CF changes with h1. I finally solved that. And with Newton's
method (also known as the Newton–Raphson method, the solution is quick and accurate.
Can anyone show the derivative or show a better method?
posted by Ragknot at
6:11 AM
12 Comments:
The derivative needed for h1
Given:
dh=3
bw=12
cr=4
cs=0.01
ss1=6
ss2=6
ss3=6
CF=1.5
ac=((bw^2*cs+2*bw*(dh-h1)*(1+cs*ss2)+(dh-h1)^2*(1+cs*ss2)*(ss2+ss3))/(2-2*cs*ss3))
ad=((cr^2*cs+h1^2*(ss1+ss2)*(1+cs*ss2)+2*cr*(h1+cs*h1*ss2))/(2-2*cs*ss1))
CF=ac/ad
ac=A1*h1^2+B1*h1+C1
(dh-h1)^2=dh^2-2*dh*h1+h1^2
A1 = (1 + cs * ss2) * (ss2 + ss3) / (2 - 2 * cs * ss3)
B1 = (-2 * bw * (1 + cs * ss2) - 2 * dh * (1 + cs * ss2) * (ss2 + ss3)) / (2 - 2 * cs * ss3)
C1 = (dh ^ 2 * (1 + cs * ss2) * (ss2 + ss3) + dh * 2 * bw * (1 + cs * ss2) + (bw ^ 2 * cs)) / (2 - 2 * cs * ss3)
ad=A2*h1^2+B2*h1+C2
A2 = (ss1 + ss2) * (1 + cs * ss2) / (2 - 2 * cs * ss1)
B2 = 2 * cr * (1 + cs * ss2) / (2 - 2 * cs * ss1)
C2 = cr ^ 2 * cs / (2 - 2 * cs * ss1)
(A1-CF*A2)*h^2+(B1-CF*B2)*h1+(C1-CF*C2)=0
a = A1 - CF * A2
b = B1 - CF * B2
c = C1 - CF * C2
plug into quadratic equation for solution
r1 = (-b + (b ^ 2 - 4 * a * c) ^ 0.5) / (2 * a)
r2 = (-b - (b ^ 2 - 4 * a * c) ^ 0.5) / (2 * a)
in this case r1=-19.54460, r2=1.54460
presumably -ve solutions are to be rejected and if b^2-4*ac <0 the solutions are imaginary and no real solution exists.
r2=1.54460 is the desired solution.
No derivative or iterative root finding techniques required. Quadratic equation gives the value of h1.
Recap:
A1 = (1 + cs * ss2) * (ss2 + ss3) / (2 - 2 * cs * ss3)
B1 = (-2 * bw * (1 + cs * ss2) - 2 * dh * (1 + cs * ss2) * (ss2 + ss3)) / (2 - 2 * cs * ss3)
C1 = (dh ^ 2 * (1 + cs * ss2) * (ss2 + ss3) + dh * 2 * bw * (1 + cs * ss2) + (bw ^ 2 * cs)) / (2 - 2 * cs * ss3)
A2 = (ss1 + ss2) * (1 + cs * ss2) / (2 - 2 * cs * ss1)
B2 = 2 * cr * (1 + cs * ss2) / (2 - 2 * cs * ss1)
C2 = cr ^ 2 * cs / (2 - 2 * cs * ss1)
a = A1 - CF * A2
b = B1 - CF * B2
c = C1 - CF * C2
r1 = (-b + (b ^ 2 - 4 * a * c) ^ 0.5) / (2 * a)
r2 = (-b - (b ^ 2 - 4 * a * c) ^ 0.5) / (2 * a)
Cam
That looks very good.
The "1.54460" is close.
Testing it gives a
CF = 1.500013139, which tells me it right.
I am going examine your methods.
The derivative I worked out is about 4 lines long... I thought there should be an easier way.
Usually the response is "You got it right." But I have to say you beat the "right" answer with a better one.
Cam,
I had 4 methods to evaluate.
Each trial is for 100,000 solutions
Your method is case 5.
You beat the heck out of me.
Congrates
'Case: 1 Last h1: 1.54460427163025 Second: 1.71875
'Case: 2 Last h1: 1.54460427163025 Second: 6.195313
'Case: 3 Last h1: 1.54460427163025 Second: 11.52344
'Case: 4 Last h1: 1.54460427163025 Second: 13.35156
'Case: 5 Last h1: 1.54460427163026 Second: 0.953125
If you still want the derivative, here it is:
((2-2*cs*ss1)*(-2*(dh-h1)*(cs*ss2+1)*(ss3+ss2)-2*bw*(cs*ss2+1)))/((2*cr*(cs*h1*ss2+h1)+h1^2*(ss2+ss1)*(cs*ss2+1)+cr^2*cs)*(2-2*cs*ss3))-((2-2*cs*ss1)*
(2*h1*(ss2+ss1)*(cs*ss2+1)+2*cr*(cs*ss2+1))*
((dh-h1)^2*(cs*ss2+1)*(ss3+ss2)+2*bw*(dh-h1)*(cs*ss2+1)+bw^2*cs))/(
(2*cr*(cs*h1*ss2+h1)+h1^2*(ss2+ss1)*(cs*ss2+1)+cr^2*cs)^2*(2-2*cs*ss3))
Right Cam,
I checked and it works, I'll add it to my time test soon.
You will be case 6
Case: 1 Last h1: 1.63367113963684 Second: 0.734375
Case: 2 Last h1: 1.63367113963684 Second: 0.625
Case: 3 Last h1: 1.63367113963684 Second: 1.0390625
Case: 4 Last h1: 1.63367113963684 Second: 1.359375
Case: 5 Last h1: 1.63367113963684 Second: 0.1015625
Case: 6 Last h1: 1.63367113963684 Second: 0.6640625
Ross, good job, but I notice that my derivative is slightly different from mine. If give a different value, but it still works.
I'll give you mine and see if you see what the diff is.
oh, I mean mine differs from yours.
I see mine is twice as long. I think you did a good job of simplifying it.
fxlg=(((h1^2-dh^2)*(1+cs*ss2)*(ss2+ss3)-dh*(2*bw)*(1+cs*ss2)-(bw^2*cs))*(2*cr)*(cs*ss2)+((h1^2-dh^2)*(1+cs*ss2)*(ss2+ss3)-dh_
*(2*bw)*(1+cs*ss2)-(bw^2*cs))*(2*cr)+((2*dh*h1^2-2*dh^2*h1)*(1+cs*ss2)*(ss2+ss3)+(h1^2-2*dh*h1)*(2*bw)*_
(1+cs*ss2)-2*h1*(bw^2*cs))*(ss1+ss2)*(1+cs*ss2)+((2*h1-2*dh)*(1+cs*ss2)*(ss2+ss3)-(2*bw)*(1+cs*ss2))*(cr^2*cs))*(2-2*cs*ss1)/(h1^2*(2-2*cs*ss3)*(2*cr)^2*_
(cs*ss2)^2+(2*h1^2*(2-2*cs*ss3)*(2*cr)^2+(2*h1^3*(2-2*cs*ss3)*(ss1+ss2)*(1+cs*ss2)+2*h1*(2-2*cs*ss3)*(cr^2*cs))*(2*cr))*(cs*ss2)+h1^2*_
(2-2*cs*ss3)*(2*cr)^2+(2*h1^3*(2-2*cs*ss3)*(ss1+ss2)*(1+cs*ss2)+2*h1*(2-2*cs*ss3)*(cr^2*cs))*(2*cr)+h1^4*(2-2*cs*ss3)*(ss1+ss2)^2*(1+cs*ss2)^2+2*_
h1^2*(2-2*cs*ss3)*(cr^2*cs)*(ss1+ss2)*(1+cs*ss2)+(2-2*cs*ss3)*(cr^2*cs)^2)
lol. I just plugged it into wxmaxima, don't remember hitting "simplify" but maybe I did.
Ragknot, you had posted the following in a different thread
"Cam,
I need to check before r1 and r2 are computed.
If Cf=1 then
.......h1 = -c / b
else
.......r1 =...
.......r2 =...etc.
If CF = 1 then A will be zero.
Thanks again."
before you compute r1 and r2 you need to check that a!=0.
in this case that occurs when CF=1, because ss1 and ss3 are equal, but if you change the ratio of ss1/ss3 the condition a=0 will occur at a different CF, resulting in a divide by 0 error.
if a=0 the equation is
0*h1^2+b*h1+c=0
or
b*h1+c=0
solving for h1....
h1=-c/b
Cam
I have...
Function camver(dh, bw, cr, cs, ss1, ss2, ss3, cf)
Dim a1, b1, c1, a2, b2, c2, a, b, c, r1, r2, r3 As Double
a1 = (1 + cs * ss2) * (ss2 + ss3) / (2 - 2 * cs * ss3)
b1 = (-2 * bw * (1 + cs * ss2) - 2 * dh * (1 + cs * ss2) * (ss2 + ss3)) / (2 - 2 * cs * ss3)
c1 = (dh ^ 2 * (1 + cs * ss2) * (ss2 + ss3) + dh * 2 * bw * (1 + cs * ss2) + (bw ^ 2 * cs)) / (2 - 2 * cs * ss3)
a2 = (ss1 + ss2) * (1 + cs * ss2) / (2 - 2 * cs * ss1)
b2 = 2 * cr * (1 + cs * ss2) / (2 - 2 * cs * ss1)
c2 = cr ^ 2 * cs / (2 - 2 * cs * ss1)
a = a1 - cf * a2
b = b1 - cf * b2
c = c1 - cf * c2
If a <> 0 Then
r1 = (-b + (b ^ 2 - 4 * a * c) ^ 0.5) / (2 * a)
r2 = (-b - (b ^ 2 - 4 * a * c) ^ 0.5) / (2 * a)
Else: r3 = -c / b
End If
If r3 > 0 Then camver = r3
If r1 > 0 Then camver = r1
If r2 > 0 Then camver = r2
End Function
Prior to calculating r1 and r2 , I would suggest checking to see if (b^2-4ac)>=0. If it is <0 the the equation only has imaginary solutions and you will be trying to calculate a -ve square root, leading to an error.
Additionally, without looking too hard at the equations, it is plausible that r1 and r2 could both be positive yielding two viable solutions (I have not confirmed this). (If b^2-4ac=0 then the two solutions would be identical) You may want to return both solutions, if both are positive, or return only one solution based on some criteria (smallest, largest,etc.)
Cam
Post a Comment
Links to this post:
Create a Link
<< Home