Unstable Simultaneous Linear Equations!
monir
1) The following set of simultaneous linear equations produces
drastically different results depending on the method and on the order
of the equations in the set!! Strange!
The set can't be described as a true sparse system!! Or, Is it ??
For ease in typing, let me use 4 equations:
ax = b
a = 1.0, -0.026625, 0.00070892, -1.8875E-05 ,b= 0.99814
a = 0.0, 1.000000, -0.05325100, 0.00212670 ,b= -0.13086
a = 1.0, -0.024859, 0.00061796, -1.5362E-05 ,b= 0.98479
a = 0.0, 1.000000, -0.04971800, 0.00185390 ,b= 132.5085
2) Here're some of the solutions, none of them is correct based on
back substitution into ALL the 4 equations.
3) One solution:(PlotIT package)
x1 = -11.911592173888
x2 = -2614.346224874
x3 = -142412.54647574
x4 = -2337438.9216864
4) Another solution: (Fortran 77 Gaussian elimination )
x1 = 0.0
x2 = -4.407E+01
x3 = 0.0
x4 = 2.338E+04
5) A third solution:(Excel VBA, Gauss method with the above equations
2 and 3 interchanged)
x1 = -1.201E-15
x2 = -2.842E-14
x3 = 4.205E+03
x4 = 1.051E+05
6) I suspect other mathematical packages (Maple, Mathematica, etc.)
would produce even different solution for the same problem!
Is there a trick to manipulate the equations in a specific way to
eliminate this type of instability ??
Thank you kindly for your expert advice.
Monir
Evgenii Rudnyi
SciPy, I obtain
>>>print svdvals(A)
[ 1.43349551e+00 1.39704245e+00 2.50232970e-03 3.40713090e-09]
This means that numerically you have a matrix with a rank of 3, that
is, 3 equations with 4 unknowns. You may want to solve this as the
least-squared problem (lstsq in SciPy).
Evgenii
http://MatrixProgramming.com
Helmut Jarausch
> Hello;
>
> 1) The following set of simultaneous linear equations produces
> drastically different results depending on the method and on the order
> of the equations in the set!! Strange!
> The set can't be described as a true sparse system!! Or, Is it ??
> For ease in typing, let me use 4 equations:
> ax = b
> a = 1.0, -0.026625, 0.00070892, -1.8875E-05 ,b= 0.99814
> a = 0.0, 1.000000, -0.05325100, 0.00212670 ,b= -0.13086
> a = 1.0, -0.024859, 0.00061796, -1.5362E-05 ,b= 0.98479
> a = 0.0, 1.000000, -0.04971800, 0.00185390 ,b= 132.5085
>
As Evgenii Rudnyi already told you, your system is ill-conditioned.
Where does this system come from? Are you trying to interpolate data
with an (unsuitable) model?
This is no "trick" to overcome these problems.
Either you come up a different system (i.e. coming from a different model)
or you have to use higher precision - but still the solution will
be very sensitive to even tiny perturbations in either the coefficients of
the matrix or the right hand side. Thus, if these come from measurements,
it's useless to use higher precision.
--
Helmut Jarausch
Lehrstuhl fuer Numerische Mathematik
RWTH - Aachen University
D 52056 Aachen, Germany
aruzinsky
> ...
> Is there a trick to manipulate the equations in a specific way to
> eliminate this type of instability ??
You can use regularized least squares (LS) which, in effect,
introduces fictitious data and a small bias to the solution. For
example, you can solve
min || Bx - y ||
where
B = [ AT c*I ]T
y = [ bT 0T ]T
c = the smallest constant for which the solution appears stable
The "smallest" c will depend on the accuracy of the LS algorithm and I
suspect that QR decomposition by givens rotation is most accurate.
JCH
"monir" <mon...@mondenet.com> schrieb im Newsbeitrag
news:2609728c-8cbc-491c...@c33g2000hsd.googlegroups.com...
I tried with a program I've written myself.
Input data:
n = 4
A*x
1
0
1
0
-0,026625
1
-0,024859
1
0,00070892
-0,05325100
0,00061796
-0,04971800
-0,000018875
0,00212670
-0,000015362
0,00185390
b
0,99814
-0,13086
0,98479
132,5085
Result:
x( 1) = 480,3657
x( 2) = 54903,42
x( 3) = 2095408
x( 4) = 2,665117E+07
Error = b_cal - b
error( 1) = 2E-26
error( 2) = 0
error( 3) = 2,9E-23
error( 4) = 0
Using 128 bit values (28 decimals)
--
Regards/Grüße http://home.arcor.de/janch/janch/menue.htm
Jan C. Hoffmann eMail aktuell: ja...@nospam.arcornews.de
Microsoft-kompatibel/optimiert für IE7+OE7
monir
> "monir" <mon...@mondenet.com> schrieb im Newsbeitragnews:2609728c-8cbc-491c...@c33g2000hsd.googlegroups.com...
>
> - Show quoted text -
Hello;
Your results:
> x( 1) = 480,3657
> x( 2) = 54903,42
> x( 3) = 2095408
> x( 4) = 2,665117E+07
do not satisfy the original equations!
I'm assuming your "," is a decimal point "." Correct ??
Monir
dave
algorithms. If I use the LU algorithm with partial pivoting (see
numerical recipes) and
double precision arithmetic I get
y= 480.36564993 54903.4222815 2095408.34306 26651166.8476
with an error |a*y -b| = 5.5825e-12
monir
>
> - Show quoted text -
Hello;
Thank you all for your replies.
Let me go back one step with the hope it might trigger some additional
thoughts!
1) I've a tiny "gap" in the analytical derivation between two points
P1 and P2.
The (x, y) coordinates and the 1st derivative are known at each of
these two point:
P1:: ( -0.026625, 0.998145) and (dy/dx) = -0.1308614
P2:: ( -0.024859, 0.984786) and (dy/dx) = 132.508494
2) The easiest and simplest approach, I thought, would be to join the
two points by a polynomial:
y = A + Bx + Cx^2 + Dx^3
with the poly. coeff. determined by satisfying the 4 conditions
described in item 1. above and solving the resulting 4 SLE.
Looks too simple to cause such difficulty! The approach (or similar!)
if successful, could easily be incorporated in the analytical model
with no much overhead!!
Any suggestions ?? Thank you.
Monir
Helmut Jarausch
a badly conditioned system.
dave
I think you are missing the point. the answer is the same as that
obtained using 128 bit arithmetic. My point is that in this case the
extra precision is not necessary if one uses a better algorithm.
JCH
Hello;
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
Yes, the program works with decimal commas. Replace data with points.
monir
>
> - Show quoted text -
Dave;
>If I use the LU algorithm with partial pivoting (see numerical recipes) and double precision
>arithmetic I get:
>y= 480.36564993 54903.4222815 2095408.34306 26651166.8476
>with an error |a*y -b| = 5.5825e-12
Unfortunately, the above LU results don't satisfy the original
equations.
Monir
dave
Yes they do unless you are multiplying incorrectly. However
one probably needs to use the full precision of the solution.
For what I submitted here the differencie is about 1.e-9.
Either that or you are doing something wrong like using the transpose of
the matrix instead of the matrix.
Evgenii Rudnyi
> >> arithmetic I get:
> >> y= 480.36564993 54903.4222815 2095408.34306 26651166.8476
> >> with an error |a*y -b| = 5.5825e-12
>
> > Unfortunately, the above LU results don't satisfy the original
> > equations.
> > Monir
>
> Yes they do unless you are multiplying incorrectly. However
> one probably needs to use the full precision of the solution.
> For what I submitted here the differencie is about 1.e-9.
> Either that or you are doing something wrong like using the transpose of
> the matrix instead of the matrix.
I have computed the norm |a*x - b| as a function of the number of
significant digits in x (n below)
n norm
5 0.624925281878
6 0.00290674506931
7 0.000365738398617
8 0.000309398615717
9 3.73063271081e-005
10 3.90003295144e-006
11 6.9922907766e-007
12 6.44133717916e-008
13 7.98163284379e-009
14 5.46330840209e-010
15 1.46064849646e-011
As you see the solution is unstable. It is hardly possible to use it
in practice. It is necessary to change the original problem. When the
ratio of singular values is about 1e10, it does not make to much sense
to proceed. As I have said, the numerical rank is 3.
dave
different packages produced different solutions. Then a poster gave the
correct solution using 128 byte arithmetic. All I did was to point out
that one could get the correct solution by using a reasonably stable
algorithm with 64 bit arithmetic. The question of whether this answer is
useful in real life is something else altogether.
Dave Dodson
> Let me go back one step with the hope it might trigger some additional
> thoughts!
>
> 1) I've a tiny "gap" in the analytical derivation between two points
> P1 and P2.
> The (x, y) coordinates and the 1st derivative are known at each of
> these two point:
>
> P1:: ( -0.026625, 0.998145) and (dy/dx) = -0.1308614
> P2:: ( -0.024859, 0.984786) and (dy/dx) = 132.508494
>
> 2) The easiest and simplest approach, I thought, would be to join the
> two points by a polynomial:
> y = A + Bx + Cx^2 + Dx^3
> with the poly. coeff. determined by satisfying the 4 conditions
> described in item 1. above and solving the resulting 4 SLE.
> Looks too simple to cause such difficulty! The approach (or similar!)
> if successful, could easily be incorporated in the analytical model
> with no much overhead!!
>
> Any suggestions ?? Thank you.
Use cubic Hermite interpolation.
Let x1, y1, and d1 be the data for P1,
and x2, y2, and d2 be the data for P2.
Given a point x, x1 <= x <= x2,
let s = (x - a) / (b - a), t = s^2, and u = s^3. Then
p(x) = y1 * (2u - 3t + 1) + y2 * (-2u + 3t)
+ d1 * (u - 2t + s) + d2 * (u - t)
Just evaluate it like this. Don't try to get A, B, C, and D so that
p(x) = Ax^3 + Bx^2 + Cx + D, because that computation will be as
unstable as your system of equations.
Dave
monir
>
> - Show quoted text -
Hello;
Thank you all for your comments and suggestions.
1) The problem is clearly associated with the ill-conditioned system.
As highlighted in the OP, and since different packages produce
different results, I personally don't believe that neither increasing
the "precision" nor using a "better" algorithm would yield the correct
results.
2) That said, here's some additional info with the hope it might
trigger some additional thoughts!
(I apologize for the repetition, but it seems that my reply earlier
today has been lost among the many responses!)
3) I've a tiny "gap" in the analytical model between two points P1 and
P2.
4) The (x, y) coordinates and the 1st derivative are known at each of
these two point:
P1:: ( -0.026625, 0.998145) and (dy/dx) = -0.1308614
P2:: ( -0.024859, 0.984786) and (dy/dx) = 132.508494
5) The easiest and simplest approach, I thought, would be to join the
two points by a polynomial:
y = A + Bx + Cx^2 + Dx^3
with the poly. coeff. determined by satisfying the 4 conditions
described in item 4. above and solving the resulting 4 SLE.
Looks too simple to cause such difficulty! The approach (or similar!)
if successful, could easily be incorporated into the analytical model
with no much overhead!!
How would you mathematically formulate the situation (3, 4, 5 above)
differently ??
Any suggestions ?? Thank you.
Monir
monir
>
> - Show quoted text -
Dave Dodson;
Your suggestion of using the cubic Hermite interpolation method sounds
very reasonable and will try it shortly.
Do you know off hand the values of "a" and "b" ??
Are they "x1" and "x2" respectively ??
Thank you.
Monir
Dave Dodson
> Your suggestion of using the cubic Hermite interpolation method sounds
> very reasonable and will try it shortly.
> Do you know off hand the values of "a" and "b" ??
> Are they "x1" and "x2" respectively ??
I'm sorry. I changed notation when typing and didn't get it all
correct. Yes. a = x1 and b = x2. I also left out a factor of (x2 - x1)
in the formula for p(x). If I haven't made another mistake, the method
should read:
Given a point x, x1 <= x <= x2,
let s = (x - a) / (x2 - x1), t = s^2, and u = s^3. Then
p(x) = y1 * (2u - 3t + 1) + y2 * (-2u + 3t)
+ [d1 * (u - 2t + s) + d2 * (u - t)] * (x2 - x1)
Dave
monir
Dave;
Having made a number of attempts with Hermite cubic interpolation, the
interpolation results between P1 and P2 suggest a typo in the formula
for p(x).
Starting at point P1, p(x) goes way beyond y2 before reversing to P2,
and the conditions d1 at P1 and d2 at P2 don't appear to be
satisfied.
Will continue checking my implementation of the method.
Regards.
Monir
Dave Dodson
> Having made a number of attempts with Hermite cubic interpolation, the
> interpolation results between P1 and P2 suggest a typo in the formula
> for p(x).
> Starting at point P1, p(x) goes way beyond y2 before reversing to P2,
> and the conditions d1 at P1 and d2 at P2 don't appear to be
> satisfied.
> Will continue checking my implementation of the method.
Monir, I e-mailed you an Excel spreadsheet showing the algorithm.
Dave
monir
Dave;
Thank you very much. I've just emailed you some comments. Will post
here the final resolution.
Monir
Helmut Jarausch
You don't need a linear at all.
Just use the so-called form functions.
In Scilab Notation
function y=N1(x)
y=(1+2*x)*(1-x)^2
endfunction
function y=N2(x)
y=x*(1-x)^2
endfunction
function y=N3(x)
y= x^2*(3-2*x)
endfunction
function y=N4(x)
y= (x-1)*x^2
endfunction
// N(0) N'(0) N(1) N'(1)
// N1 1 0 0 0
// N2 0 1 0 0
// N3 0 0 1 0
// N4 0 0 0 1
P1 =[-0.026625; 0.998145]
P1S= -0.1308614
P2=[-0.024859; 0.984786]
P2S=132.508494
x1= P1(1); y1= P1(2);
x2= P2(1); y2= P2(2);
dx=x2-x1
function y=T(x)
y= (x-x1)/dx
endfunction
function y=N1T(x)
y=N1(T(x))
endfunction
function y=N2T(x)
y=N2(T(x))*dx
endfunction
function y=N3T(x)
y=N3(T(x))
endfunction
function y=N4T(x)
y=N4(T(x))*dx
endfunction
// N(x1) N'(x1) N(x2) N'(x2)
// N1T 1 0 0 0
// N2T 0 1 0 0
// N3T 0 0 1 0
// N4T 0 0 0 1
function y=yourfunction(z)
y= y1*N1T(z) + y2*N3T(z) + P1S*N2T(z) + P2S*N4T(z)
endfunction
monir
>
> - Show quoted text -
Helmut;
Clever idea!
If you still have the work file, what y-value you get for x =
-0.0255 ??
Thanks again.
Monir
Dave Dodson
>
> - Show quoted text -
Monir, if you look closely, you will see that what Helmut proposes is
the same thing I proposed. His T(x) is the expression I used to
compute s, and his functions N1, N2, N3, and N4 are the same as the
four polynomials in my expression for p(x).
Furthermore, since the cubic polynomial that passes through (x1,y1)
with slope d1 and (x2,y2) with slope d2 is unique, this answer will be
0.954300184.
Dave
monir
>
> - Show quoted text -
Dave;
You're correct!
Here's a summary so far:
1) The coordinates (x, y) and the 1st derivative "d" at two points P1
and P2 are known.
General:
x1 < x2 < 0
y1 > y2 > 0
Sign of d1 and d2 could be any of the 4 +/- combinations.
For example:
..........x.........y.........d=(dy/dx)
P1:: -0.026625 0.998145 -0.1308614
P2:: -0.024859 0.984786 132.5084940
Goal:
To join points P1 and P2 by a smooth exact (not LS) curve with no
extremum and no inflection points in the range x1 - x2
(By plotting the two points and the slope lines, it appears a solution
is clearly feasible!)
2) The easiest and simplest approach, I thought, would be to join the
two points by a polynomial:
y = A + Bx + Cx^2 + Dx^3
with the poly. coeff. determined by satisfying the 4 conditions
described in item 1. above (2 points and two slopes) and solving the
resulting 4 SLE.
The approach (according to: OP, Evgenii Rudnyi, Hulmet Jarausch, JCH,
and dave)
fails due to the inherent numerical instability in the system no
matter how clever the algorithm is or how extra precision is applied
in solving the ill-conditioned system.
3) The 2nd approach uses Hermite cubic interpolation formula
(according to Dave Dodson and Hulmet Jarausch).
The method is based on a 3rd deg poly and introduces extremum and
inflection points between points P1 and P2, and thus does not produce
the desired interpolation curve.
4) If one re-examines the problem description in item 1. above, it is
reasonable to assume that there must be a mathematical function
flexible enough to satisfy the conditions and constraints! Correct ??
Perhaps a Cartesian spiral of some kind ?? or a non-standard truncated
parabolic function ?? or a semicubical ?? or ......
Any additional thoughts ?? Thank you.
Monir
Dave Dodson
There is not a _differentiable_ function flexible enough to satisfy
the conditions and constraints for all four combinations of
derivatives. Given any curve y = f(x), according to the Mean Value
Theorem, there is a point (x3,y3) on the curve such that y'(x3) equals
the slope of the chord between (x1,y1) and (x2,y2). Since x2 > x1 and
y2 < y1, this slope is negative, so y'(x3) is negative. If d2 = y'(x2)
is positive, then the Intermediate Value Theorem states that there is
a point x4 where y'(x4) = 0. This is an extremum. Thus, your condition
cannot be met when d2 > 0. Furthermore, if d1 > 0, then again by the
Intermediate Value Theorem, there will be extremum between x1 and x3.
If d1 > 0 and d2 > 0, then there must be two extrema and if the
function is twice differentiable, there must be an inflection point.
In summary, here are the cases:
d1 < 0, d2 < 0: solution may be possible.
d1 < 0, d2 > 0: solution impossible.
d1 > 0, d2 < 0: solution impossible.
d1 > 0, d2 > 0: solution impossible.
Since polynomials are once and twice differentiable, the only
situation where a polynomial meets your conditions is when d1 < 0 and
d2 < 0.
Finally, note that for any function that is differentiable at x2, if
d2 is positive, then the curve must be increasing on some interval [x2-
epsilon, x2], so there will be a point to the left of x2 where the
function is less than y2.
In short, there is no function that satisfies your constraints for the
given problem.
Dave
monir
>
> - Show quoted text -
Dave;
Thank you very much for your thorough reply!
Will keep trying!
Monir
JCH
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
Testing my solution before (decimal commas!)
x( 1) = 480,36564993069834577726219468
x( 2) = 54903,422281484361120288636181
x( 3) = 2095408,3430612314008249198507
x( 4) = 26651166,847636842152179043019
eps = A * x - b
eps( 1) = 0,00000000000000000000000002
eps( 2) = 0,000000000000000000000000
eps( 3) = -0,00000000000000000000000006
eps( 4) = -0,000000000000000000001106
and comparing after small perturbations of b values
x_test( 1) = 480,84601558062904412303945690
x_test( 2) = 54958,325703765845481408924813
x_test( 3) = 2097503,7514042926322257447704
x_test( 4) = 26677818,014484478994331222061
Stability Test: b_test = 1,001 * b
eps = A * x_test - b_test
eps( 1) = 0,00000000000000000000000000
eps( 2) = 0,000000000000000000000000
eps( 3) = 0,00000000000000000000000000
eps( 4) = -0,000000000000000000001098
shows stability.
See also page 33
http://www.matrixanalysis.com/Chapter1.pdf
monir
> "monir" <mon...@mondenet.com> schrieb im Newsbeitragnews:3e3925c1-5008-4255...@i12g2000prf.googlegroups.com...
>
> --
> Regards/Grüße http://home.arcor.de/janch/janch/menue.htm
> Jan C. Hoffmann eMail aktuell: ja...@nospam.arcornews.de
> Microsoft-kompatibel/optimiert für IE7+OE7
Hi Jan;
Thank you once again for your thorough analysis of the problem.
1) For the same given example:
..........x.........y........d=(dy/dx)
P1:: -0.026625 0.998145 -0.1308614
P2:: -0.024859 0.984786 132.5084940
ax = b
a = 1.0, -0.026625, 0.00070892, -1.8875E-05 ,b= 0.99814
a = 0.0, 1.000000, -0.05325100, 0.00212670 ,b= -0.13086
a = 1.0, -0.024859, 0.00061796, -1.5362E-05 ,b= 0.98479
a = 0.0, 1.000000, -0.04971800, 0.00185390 ,b= 132.5085
2) I applied your solution in Excel worksheet (max allowed is double
precision):
x( 1) = 480.365 649 930 698 345 777 262 194 68
x( 2) = 54,903.422 281 484 361 120 288 636 181
x( 3) = 2,095,408.343 061 231 400 824 919 850 7
x( 4) = 26,651,166.847 636 842 152 179 043 019
3) Here's the a comparison between the input values and your results
at points P1 and P2:
......|.......input......|.....JCH results...|...input.....|....JCH
results...|
...........x........y...............y...........d=(dy/dx)........(dy/
dx)......
P1:: -0.026625 0.998145 0.957564182689907 -0.1308614
1.24646690826194
P2:: -0.024859 0.984786 1.002504622037070 132.5084940
132.81214634587400
Considerable discrepancy in the values of "y" and 1st derivative !!
4) I've also plotted the interpolation curve between P1 and P2 based
on your results. The curve is substantially off target.
5) What is encouraging, however, is that I've just and by accident
developed a simple procedure which appears to be working fine and as
desired. My elliptic interpolation procedure:
a. is numerically stable and does not require extra precision;
a. joins the two points P1 and P2 by an exact smooth continuous well-
defined curve;
b. satisfies the 4 conditions (2 points and 2 slopes);
c. produces correct interpolated y-value between y2 and y1 over the
range x2 to x1; and
d. presents no extremum or inflection points in the x-range-of-
interest.
I'm checking it carefully before posting.
Regards.
Monir
The Phantom
wrote:
Consider the first row of your matrix a. The 1,2 element is -.026625 and
the 1,3 element is shown as .00070892. But the 1,3 element should be the
square of the 1,2 element, namely .000708890625. The 1,4 element should be
the cube of the 1,2 element, namely -1.88742128906E-5. Similarly for the
elements of the 3rd row. You need to carry more digits in the derived
elements of the matrix.
You also show a few more digits for the elements of your b (column) matrix
just above. You should incorporate these in your b matrix.
If you do all this, then the coefficients for your 3rd degree polynomial
are:
A = 846.799620372
B = 97712.4903969
C = 3760588.91028
D = 48216720.1782
Evaluate your polynomial at the endpoints using these coefficients and you
will get much better results.
monir
>
> - Show quoted text -- Hide quoted text -
>
> - Show quoted text -
Hello;
1) Given:
..........x.........y........d=(dy/dx)
P1:: -0.026625 0.998145 -0.1308614
P2:: -0.024859 0.984786 132.5084940
2) The Hermite cubic interpolation formula introduces extremum and
inflection points between points P1 and P2, and thus does not produce
the desired interpolation curve.
3) The JCH SLE solution (applied in Excel w/s, double precision):
y = A + B.x +C.x^2 + D.x^3
A = 480.365 649 930 698 345 777 262 194 68
B = 54 903.422 281 484 361 120 288 636 181
C = 2 095 408.343 061 231 400 824 919 850 7
D = 26 651 166.847 636 842 152 179 043 019
Comparison between the input values and JCH results (truncated) at the
endpoints P1 and P2:
....|........input......|..JCH ...|...input....|..JCH results
..........x.........y.........y......(dy/dx).......(dy/dx)
P1:: -0.026625 0.998145 0.957564 -0.1308614 1.2464669
P2:: -0.024859 0.984786 1.002504 132.5084940 132.8121463
Considerable discrepancy in the values of "y" and 1st derivative!
4) Your (The Phantom) SLE solution (which're drastically different
from JCH's):
A = 846.799620372
B = 97 712.4903969
C = 3 760 588.91028
D = 48 216 720.1782
and at the endpoints:
y1 = 0.99814391 ... (dy/dx) = 2.35924858
y2 = 0.98481476 ... (dy/dx) = 133.04470474
Still considerable discrepancies in the value of 1st derivative at
P1!!
So it appears that extra precision wouldn't solve the problem.
5) Here's a possible solution.
If we join the two points P1 and P2
...........x.........y........d=(dy/dx)
P1:: -0.0266255 0.9981449 -0.1593857
P2:: -0.0248588 0.9847863 132.5084940
by the elliptic equation:
Ax^2 + By^2 + Cx + Dy = 1
and solve the resulting 4 SLE:
0.00070892*A + 0.99629333*B - 0.0266255*C + 0.99814495*D = 1
0.00061796*A + 0.96980405*B - 0.0248588*C + 0.98478630*D = 1
-0.05325094*A + 0.31817996*B + 1.0000000*C - 0.15938565*D = 0.0
-0.04971757*A + 260.985098*B + 1.0000000*C - 132.508494*D = 0.0
we get:
A = - 48.857138
B = - 0.994123
C = -2.605710
D = 1.959330
6) The results at the endpoints are:
....|.........input.........|.elliptic.|...input...|..elliptic
...........x..........y..........y.......(dy/dx).....(dy/dx)
P1:: -0.0266255 0.9981449 0.9981449 -0.1593857 -0.1593857
P2:: -0.0248588 0.9847863 0.9847863 132.5084940 132.508494
7) The full range Elliptic interpolation results for x1 <= x <= x2:
X...........Y results...(dY/dX)results
-0.0266255 0.9981449 -0.1593857
-0.0265812 0.9981341
-0.0265370 0.9981157
-0.0264928 0.9980896
-0.0264486 0.9980558
-0.0264043 0.9980143
-0.0263601 0.9979650
-0.0263159 0.9979077
-0.0262716 0.9978425
-0.0262274 0.9977690
-0.0261832 0.9976873
-0.0261389 0.9975971
-0.0260947 0.9974983
-0.0260505 0.9973906
-0.0260062 0.9972737
-0.0259620 0.9971475
-0.0259178 0.9970116
-0.0258736 0.9968656
-0.0258293 0.9967093
-0.0257851 0.9965420
-0.0257409 0.9963634
-0.0256966 0.9961728
-0.0256524 0.9959696
-0.0256082 0.9957531
-0.0255639 0.9955224
-0.0255197 0.9952765
-0.0254755 0.9950142
-0.0254313 0.9947341
-0.0253870 0.9944346
-0.0253428 0.9941137
-0.0252986 0.9937689
-0.0252543 0.9933970
-0.0252101 0.9929941
-0.0251659 0.9925549
-0.0251216 0.9920721
-0.0250774 0.9915354
-0.0250332 0.9909289
-0.0249889 0.9902261
-0.0249447 0.9893755
-0.0249005 0.9882451
-0.0248563 0.9854568
-0.0248565 0.9852447
-0.0248568 0.9851569
-0.0248573 0.9850326
-0.0248578 0.9849374
-0.0248583 0.9848570
-0.0248588 0.9847863 132.508494
7) The above elliptic interpolation procedure works fine and as
desired:
a. is numerically stable and does not require extra precision;
b. satisfies the 4 conditions (2 points and 2 slopes);
c. accommodates any possible +/- sign combinations of d1 and d2;
d. joins the two points P1 and P2 by an exact smooth continuous well-
defined curve;
e. produces interpolated y-value bounded by y2 and y1 over the range
x2 to x1; and
f. presents no extremum or inflection points in the x-range-of-
interest.
Any comments or thoughts ?? Thank you.
Monir
Dave Dodson
The curve is not a function. Look at the last 8 lines in your table,
and you will see that x first increases and then decreases. What value
of y will you produce for x = -0.0248575, for example? Will you choose
the upper value, for which y will be between 0.9852447 and 0.9854568,
or the lower value, for which y will be between 0.9849374 and
0.9850326? These correspond to whether you choose the + sign or the -
sign in the quadratic formula when you are solving the equation of the
ellipse for y for a given value of x.
When you combine this curve with the curve to the right, you actually
will have 3 possible values for y(-0.0248575), the previous two, plus
the value from the curve to the right.
There is a vertical tangent somewhere near x = -0.0248563. At this
point, the derivative y' is undefined. Furthermore, there is a cusp
where the elipse joins the curve to the right, at (-0.0248588,
0.9847863). As you traverse the ellipse toward this point, there will
be a sharp point where you will have to reverse directions. Does these
conditions really satisfy your definition of _smooth_?
Dave
monir
>
> - Show quoted text -
Dave;
1) Circles, ellipses, parabolas, spirals, etc., are all multi-value,
exact, well-defined, smooth plane curves, with vertical tangent at
some point(s) along the curve before reversing directions, and "x"
increases and decreases as one traverses the curve.
Such chacteristics of the named smooth curves present no problems!
2) One uses +ve or -ve sqrt for "y" depending on the value of the
independent variable "x" in relation to the centre of the curve and
the coordinates system.
For x= -0.0248575 (very close to x2), the procedure calculates y =
0.98592732 just before reaching the vertical tangent to the right, and
y = 0.98498626 just after that and moving towards the endpoint P2.
Both values are correct and are on the correct smooth elliptic
interpolation curve.
3) I haven't observed a "cusp" or a "sharp point" anywhere and
certainly not at or near P2 for any of the tested cases so far! It
would be a serious setback if your observation could be substantiated.
I've enlarged the corresponding plot and zoomed on P2. There's no
"cusp" or any abnormality to be identified anywhere!
Did you observe it visually or from the analysis of the provided
numerical table ??
Could you please elaborate?
Thank you.
Monir
Dave Dodson
>2) One uses +ve or -ve sqrt for "y" depending on the value of the
>independent variable "x" in relation to the centre of the curve and
>the coordinates system.
So I give you an x. How do you know which sign, and therefore which y
to produce?
> 3) I haven't observed a "cusp" or a "sharp point" anywhere and
> certainly not at or near P2 for any of the tested cases so far! It
> would be a serious setback if your observation could be substantiated.
> I've enlarged the corresponding plot and zoomed on P2. There's no
> "cusp" or any abnormality to be identified anywhere!
> Did you observe it visually or from the analysis of the provided
> numerical table ??
> Could you please elaborate?
In an early post in this thread, you stated that you had a tiny gap in
an analytic development between two points, and that you know the
points and the slopes at each. To me that implied that you have a
curve on each side of the gap, and you want to fill the gap with an
interpolating curve. If this is not the case, then my comment is
irrelevant.
Dave
monir
Dave;
Thanks again for your prompt reply and for your interest in the
problem.
> So I give you an x. How do you know which sign, and therefore which y to produce?
As I suggested earlier, there're 4 possible scenarios that one must
check after solving the 4 SLE for A, B, C, D in order to channel the
procedure to construct the correct elliptic interpolation curve:
given: point P1(x1,y1) & slope d1
point P2(x2,y2) & slope d2
scen1: IF y1 >= k & y2 >= k i.e.; both P1 and P2 on the upper surface
..........use "+" sqrt() for x1 <= x <= x2
scen2: IF y1 > k and y2 < k i.e.; P1 on the upper and P2 on the lower
surface
..........use "+" sqrt() for x1 <= x <= h+a on the upper surface
..........use "-" sqrt() for h+a > x >= x2 on the lower surface
scen3: IF y1 < k and y2 > k i.e.; P1 on the lower and P2 on the upper
surface
..........use "-" sqrt() for x1 >= x <= h-a on the lower surface
..........use "+" sqrt()for h-a < x <= x2 on the upper surface
scen4: IF y1 < k and y2 < k i.e.; P1 and P2 on the lower surface
..........use "-" sqrt() for x1 >= x <= h-a on the lower surface
..........use "+" sqrt() for h-a < x <= h+a on the upper surface
..........use "-" sqrt() for h+a > x > x2 on the lower surface
(scen2 above applies to the worked example and the solution provided
earlier)
where:
elliptic interpolation curve: Ax^2 + By^2 + Cx + Dy = 1
A = (b^2/a^2)/[b^2 - k^2 - (b^2/a^2)h^2]
B = 1/[b^2 - k^2 - (b^2/a^2)h^2]
C = -2h(b^2/a^2)/[b^2 - k^2 - (b^2/a^2)h^2]
D = -2k/[b^2 - k^2 - (b^2/a^2)h^2]
with centre at (h,k)
ellipse axis along the X-axis = 2a
ellipse axis along the Y-axis = 2b
>In an early post in this thread, you stated that you had a tiny gap in an
>analytic development between two points, and that you know the points and
>the slopes at each. To me that implied that you have a curve on each side
>of the gap, and you want >to fill the gap with an interpolating curve.
>If this is not the case, then my comment is irrelevant.
This is exactly the case!
HTH.
Monir
Dave Dodson
Here is where I see the problem with your curve not being single
valued. If I give you an x > x2, where the curve has two values, how
do you know whether I want the upper one or the lower one so that you
know whether to use the + sign or the - sign?
Okay. If you trace along the ellipse from (x1,y1), you get to the
point (x2,y2), headed downward and to the left. The curve that you are
matching function value and slope is heading upward and to the right.
This creates a cusp. If, for example, you were marching along the
ellipse from (x1,y1), you would have to do an about-face at point
(x2,y2) in order to continue marching along the analytically developed
curve that begins at that point. I wouldn't consider this "smooth." Do
you?
Now, for an x > x2, it is possible that a vertical line cuts your
combined curve in 3 places: the upper curve of the ellipse, the lower
curve of the ellipse, and the analytically-developed curve starting at
(x2,y2). How do you know which of the three possible y values I want?
Am I beating a dead horse?
Dave
monir
>
>
> - Show quoted text -
Dave;
There's a good chance that I overlooked little details in describing
the procedure in my post which's causing such confusion! This is
likely the case since I haven't had any problem (so far) in applying
the procedure for 10s of different scenarios.
Since the proof of the pudding in the eating, let me suggest the
following.
Elliptic Equation:
Ax^2 + By^2 + Cx + Dy = 1
1) Plot the numerical table provided earlier (Mar 31, 3:36 local) for:
...........x.........y........d=(dy/dx)
P1:: -0.0266255 0.9981449 -0.1593857
P2:: -0.0248588 0.9847863 132.5084940
A = -48.857138
B = -0.994123
C = -2.605710
D = 1.959330
2) Plot the following 2nd example:
..........x.........y......d=(dy/dx)
P1:: -0.057785 0.996125 0.496957
P2:: -0.048650 0.985900 43.111960
A = -2.284E+00
B = -1.021E+00
C = -2.540E-01
D = 2.014E+00
The full-range Elliptic interpolation results:
X..........Y results...(dY/dX)results
-0.057785 0.996125 0.496957
-0.057557 0.996232
-0.057328 0.996327
-0.057100 0.996409
-0.056871 0.996478
-0.056643 0.996536
-0.056414 0.996583
-0.056186 0.996618
-0.055957 0.996641
-0.055729 0.996654
-0.055500 0.996655
-0.055272 0.996644
-0.055043 0.996623
-0.054815 0.996590
-0.054586 0.996545
-0.054358 0.996490
-0.054129 0.996422
-0.053901 0.996342
-0.053672 0.996250
-0.053444 0.996145
-0.053215 0.996027
-0.052987 0.995896
-0.052758 0.995750
-0.052530 0.995590
-0.052301 0.995414
-0.052073 0.995222
-0.051844 0.995012
-0.051616 0.994783
-0.051387 0.994534
-0.051159 0.994262
-0.050930 0.993966
-0.050702 0.993642
-0.050474 0.993287
-0.050245 0.992895
-0.050017 0.992460
-0.049788 0.991971
-0.049560 0.991414
-0.049331 0.990763
-0.049103 0.989968
-0.048874 0.988904
-0.048646 0.986261
-0.048646 0.986261
-0.048646 0.986147
-0.048646 0.986100
-0.048647 0.986063
-0.048647 0.986033
-0.048648 0.986006
-0.048648 0.985982
-0.048649 0.985959
-0.048649 0.985938
-0.048649 0.985919
-0.048650 0.985900 43.111960
3) In the above two examples:
Do you actually observe any cusps, discontinuities, cross-overs, sharp
points, etc., anywhere along the "smooth" elliptic curve joining P1
and P2 ??
I don't !!
4) If your answer is YES, then there's a problem with the procedure.
If your answer is NO, then it would be clear that I overlooked some
details in describing the method and I'll try to provide a complete
description.
Regards.
Monir
Dave Dodson
As I have been trying to explain, the cusp isn't within the ellipse.
It it where the ellipse joins the curve to the right at point P2. As
you traverse the ellipse from P1 to P2, you go down and to the right
until you get to the vertical tangent somewhere beyond x2. Then you
continue down, but go to the left until you get to P2. At P2, you
switch from the ellipse to the analytically derived curve that goes up
and to the right. The cusp is where you switch from going down and to
the left to going up and to the right, exactly at P2. Oh for a
picture!
monir
>
> - Show quoted text -- Hide quoted text -
>
> - Show quoted text -
Dave;
Thank you for emailing me the file. It's true that a picture is worth
a 1,000 words!
The "cusp" you referred to on your plot has nothing to do with the
elliptic interpolation procedure. It is not a "cusp" but rather a
plotting issue!
In the provided worked examples, the selected equi-interval x values
on the upper surface near P2 is ~ 200 times that on the lower surface,
and thus the Excel "smooth line" would be hard-pressed to join the
points smoothly.
Choose x-interval on the upper and lower surfaces of a similar order
of magnitude.
(I've just emailed you an Excel file showing how the correct plot
should look like.)
Still for a better display, one could apply a cosine distribution for
the x values, resulting in a concentration of interpolated points near
the ends.
Another possibility is to use a general ellipse at an angle, but one
additional condition would be needed for an exact solution of the
problem.
Regards.
Monir
Dave Dodson
> Thank you for emailing me the file. It's true that a picture is worth
> a 1,000 words!
>
> The "cusp" you referred to on your plot has nothing to do with the
> elliptic interpolation procedure. It is not a "cusp" but rather a
> plotting issue!
I think you just don't get it. I know what a cusp is, and the graph
you sent me exhibits a cusp at the point P2 just as much as the graph
I sent you, but if you are happy with the fit, then fine.
> In the provided worked examples, the selected equi-interval x values
> on the upper surface near P2 is ~ 200 times that on the lower surface,
> and thus the Excel "smooth line" would be hard-pressed to join the
> points smoothly.
>
> Choose x-interval on the upper and lower surfaces of a similar order
> of magnitude.
> (I've just emailed you an Excel file showing how the correct plot
> should look like.)
The point I was making was about the "point" formed where the black
curve and the yellow curve join, not about the spacing of plotting
points along the black curve or how that affects the appearance of the
graph. If you think the black and yellow curves join smoothly, then
fine.
> Still for a better display, one could apply a cosine distribution for
> the x values, resulting in a concentration of interpolated points near
> the ends.
>
> Another possibility is to use a general ellipse at an angle, but one
> additional condition would be needed for an exact solution of the
> problem.
I hope that I have been of some help, but I don't think I can
contribute any more.
Dave
monir
The elliptic interpolation procedure works fine for me and exactly as
desired.
Thank you all for your help. Special thanks to Dave Dodson.
Regards.
Monir
The Phantom
<SNIP>
>
>4) Your (The Phantom) SLE solution (which're drastically different
>from JCH's):
>A = 846.799620372
>B = 97 712.4903969
>C = 3 760 588.91028
>D = 48 216 720.1782
>and at the endpoints:
>y1 = 0.99814391 ... (dy/dx) = 2.35924858
>y2 = 0.98481476 ... (dy/dx) = 133.04470474
>Still considerable discrepancies in the value of 1st derivative at
>P1!!
>So it appears that extra precision wouldn't solve the problem.
If I substitute -.026625 into the polynomial using the coefficients I gave you,
I get .998144999955 rather than the .99814391 you get. You are not using enough
precision in your arithmetic, or you didn't really use -.02662500000 as the
value of x.
I can get the value you got, if I use -.02662547000 for the x value, and using
high precision arithmetic.
The derivatives are off quite a bit because I didn't check those elements of the
matrix for correctness before I solved it.
So, I went back and redid ALL the elements of the matrix with high precision
arithmetic. Then I get for a solution:
A = 831.487177808938
B = 95 912.4414665007
C = 3 690 085.37174486
D = 47 296 658.6429856
The reason several people (including myself) got different results is because we
accepted the low precision elements you gave for your matrix. If you will
evaluate your cubic polynomial with the ABCD coefficients just above, you will
find that the y values and derivatives will be correct to 9 significant digits.
But you must use the values -.0266250000 and -.0248590000 with high precision
arithmetic.
It is possible that the values -.026625 and -.024859 are truncated or rounded
values themselves. I say that because some of your numbers were given with
differing numbers of significant digits in various places. But, if you take
those x values to be exact as shown, then the coefficients I just gave will give
you 9 digit accuracy, as I said.
The Phantom
<SNIP>
>
>5) Here's a possible solution.
>If we join the two points P1 and P2
>...........x.........y........d=(dy/dx)
>P1:: -0.0266255 0.9981449 -0.1593857
I find this confusing. In your very post, the value of dy/dx for the first
point was given as -0.13086, but now it has become -0.1593857. Which is it?
monir
>
> - Show quoted text -
Rodger;
Okay.
Given: (exact 7 significant digits)
..........x.........y........d=(dy/dx)
P1:: -0.0266255 0.9981449 -0.1308565
P2:: -0.0248588 0.9847863 132.5335679
what poly coeff A, B, C, D do you get for:
y = A + B*x +C*x^2 + D*x^3
Upon receipt, I'll provide you with an Excel plot showing the
comparison and how to construct the correct elliptic interpolation
curve.
Regards.
Monir
The Phantom
<SNIP>
>Rodger;
>
>Okay.
>Given: (exact 7 significant digits)
>..........x.........y........d=(dy/dx)
>P1:: -0.0266255 0.9981449 -0.1308565
>P2:: -0.0248588 0.9847863 132.5335679
>
>what poly coeff A, B, C, D do you get for:
>y = A + B*x +C*x^2 + D*x^3
>
>Upon receipt, I'll provide you with an Excel plot showing the
>comparison and how to construct the correct elliptic interpolation
>curve.
>
>Regards.
>Monir
I calculated the coefficients in Mathematica with 50 digit arithmetic. Other
high precision arithmetic packages could be used, such as Matlab, Maple, Scilab,
etc. Here is a link to an image showing the whole thing, including a plot of
the cubic passing through the endpoints plus the elliptic function you gave in
another post:
http://img2.freeimagehosting.net/image.php?8ada8a8e27.gif
The coefficients I got are:
A = 830.9589672389913
B = 95 850.51921961399
C = 3 687 666.010250186
D = 47 265 155.70059824
I think Excel won't be able to get this result even using double precision
arithmetic, and probably won't even be able to accurately evaluate the cubic and
the derivatives with double precision.
The results and calculations shown in the Mathematica result are quite accurate,
and can be relied on to be correct, for the starting values:
Given: (exact 7 significant digits)
..........x.........y........d=(dy/dx)
P1:: -0.0266255 0.9981449 -0.1308565
P2:: -0.0248588 0.9847863 132.5335679
This all shows that extreme care and high precision arithmetic is required to
get a good result for the cubic because the condition number for the system is
very high. The cubic is not a good function for the interpolation, but it can
be accurately calculated with high precision arithmetic.
monir
>
> - Show quoted text -
Rodger;
Thank you for the info. and for your effort.
1) I've applied your latest poly. coeffs. in Excel for comparison
purposes (file HERMITE-JCH-Dave-Rodger-Monir-Interpol-6.xls, just
emailed to you).
2) It might surprise you that your cubic poly calculations using DP
arithmetic agree extremely well with the Hermite method.
The difficulty, however and as demonstrated earlier with general cubic
poly interpolation, is the unavoidable presence of extremum and
inflection points within the x-range-of-interest, which is clearly
undesirable!
3) On the other hand, the elliptic interpolation curve appears to
satisfy all the conditions and requirements. It's critical, however,
that one applies the correct scenario in constructing the applicable
curve segment of the ellipse. The 4 possible scenarios were described
earlier (my post April 1, 12:20 am local).
The applicable scenario for the provided example is scen2
4) As a footnote, Dave Dodson has insisted that in the elliptic
interpolation procedure a "cusp" exists where the ellipse joins the
(analytical) curve at point P2.
Dave's observation suggests that although the elliptic curve and the
analytical curve share the same coordinates (x2,y2) and have the same
slope (dy/dx) at P2, the elliptic curve nevertheless has a cusp at
that point!
I haven't personally observed or been able to locate such abnormality
at or near P2 in any of the 10s test cases.
Have you observed a cusp at P2 ??
Any comments ? Thank you.
Monir
The Phantom
<SNIP>
>
>3) On the other hand, the elliptic interpolation curve appears to
>satisfy all the conditions and requirements. It's critical, however,
>that one applies the correct scenario in constructing the applicable
>curve segment of the ellipse. The 4 possible scenarios were described
>earlier (my post April 1, 12:20 am local).
>The applicable scenario for the provided example is scen2
>
>4) As a footnote, Dave Dodson has insisted that in the elliptic
>interpolation procedure a "cusp" exists where the ellipse joins the
>(analytical) curve at point P2.
>Dave's observation suggests that although the elliptic curve and the
>analytical curve share the same coordinates (x2,y2) and have the same
>slope (dy/dx) at P2, the elliptic curve nevertheless has a cusp at
>that point!
>I haven't personally observed or been able to locate such abnormality
>at or near P2 in any of the 10s test cases.
>Have you observed a cusp at P2 ??
>
>Any comments ? Thank you.
>Monir
I think what Dave is talking about is something that I also wondered about.
Your elliptic curve does have the right slope at P2, but if you follow the path
along the elliptic curve from P1 to P2, when you arrive at P2 you are going down
and to the left. Can you show some of the points in your analytical function
that are outside the interval P1 < x < P1? I assumed that as you approach P2
from the left, you would want to be moving up and to the right.
As Dave pointed out in another post, if that is the case, then if you want a
slope of 132.533+ at P2 you will have to be approaching from below and the
interpolating curve will have to have a minimum as the cubic does.
Your elliptic curls around as it approaches P2; surely that isn't what you want,
is it?
monir
>
> - Show quoted text -
The Phantom;
>"I think what Dave is talking about is something that I also wondered about."
I've recalculated the elliptic function around the vertical tangent to
the right and leading to point P2, with x-interval ~ 2.E-7
Given: (exact 7 significant digits)
..........x.........y........d=(dy/dx)
P1:: -0.0266255 0.9981449 -0.1308565
P2:: -0.0248588 0.9847863 132.5335679
pnt.......X.......Y results...(dY/dX)results
1 -0.0248585 0.9860975
2 -0.0248583 0.9860610
3 -0.0248580 0.9860221
4 -0.0248578 0.9859804
5 -0.0248575 0.9859350
6 -0.0248573 0.9858848
7 -0.0248570 0.9858278
8 -0.0248568 0.9857601
9 -0.0248565 0.9856720
10 -0.0248564 0.9856096
11 -0.0248563 0.9854591 vertical tangent
12 -0.0248563 0.9854591 vertical tangent
13 -0.0248564 0.9853086
14 -0.0248565 0.9852463
15 -0.0248568 0.9851582
16 -0.0248570 0.9850905
17 -0.0248573 0.9850335
18 -0.0248575 0.9849833
19 -0.0248578 0.9849379
20 -0.0248580 0.9848961
21 -0.0248583 0.9848573
22 -0.0248585 0.9848208
23 -0.0248588 0.9847863 132.5335679
Please identify which point has a "cusp" ?? You may refer to pnt #
instead of coordinates.
>"... and the interpolating curve will have to have a "minimum" as the cubic does."
Please identify where is the "minimum" ??
So far, I've tested the elliptic procedure for about 60 different
pairs, with different scenarios scen1, scen2, scen3, scen4, and with
different signs and magnitudes of given slopes d1 and d2, and the
procedure appears to be working fine!
The Phantom
>On Apr 3, 11:09 pm, The Phantom <phan...@aol.com> wrote:
>> On Thu, 3 Apr 2008 18:14:51 -0700 (PDT), monir <mon...@mondenet.com> wrote:
>>
>> <SNIP>
>Can you show some of the points in your analytical function
>> that are outside the interval P1 < x < P1?
You didn't answer this question and I think there may be some mis-communication
about what you want because of this very thing.
>I assumed that as you approach P2
>> from the left, you would want to be moving up and to the right.
>>
>
These points are on the very right end of the ellipse and there is no cusp
associated with this curve BY ITSELF.
>
>>"... and the interpolating curve will have to have a "minimum" as the cubic does."
>Please identify where is the "minimum" ??
The top half of the elliptic curve which connects P1 to P2 doesn't have a
minimum of the sort I was talking about.
>
>So far, I've tested the elliptic procedure for about 60 different
>pairs, with different scenarios scen1, scen2, scen3, scen4, and with
>different signs and magnitudes of given slopes d1 and d2, and the
>procedure appears to be working fine!
>
>Any comments or thoughts ?? Thank you.
>Monir
I have assumed all along that there more points, or another "analytical" curve,
to the left of P1 and to the right of P2, and that the curve you are attempting
to derive should "fill in the gap", meeting the points (or curve) outside the
interval P1 < x < P1, matching in position and slope at P1 and P2.
Imagine that you start at P1 and move on some curve toward P2. When you arrive
at P2 I assume that you want to be moving on a slope of ~132, upward and to the
right, NOT downward and to the left.
You haven't said which one of these you want, but I assumed you want to be
heading northeast, not southwest.
If you use the ellipse as you have shown, when you reach P2, you will be heading
southwest; the slope will be ~132, but you will be heading in a direction that
takes you back in the the x-interval BETWEEN P1 and P2, rather than out to the
right. You will have reversed your travel in the x-coordinate.
As I say, I have assumed that you want to be going northeast when you reach P2.
If you DO want to be going northeast, then if you travel from P1 to P2 via the
elliptic curve, you must reverse direction when you reach P2 in order to
continue out of the interval in a northeasterly direction.
Have a look at this graphic:
http://img2.freeimagehosting.net/image.php?13225904c8.gif
The top graph shows the red elliptic curve going from P1 to P2, but when you get
to P2, if you want to continue to leave the P1-P2 interval in a northearterly
direction, you must reverse direction and follow the black line, which is just a
line with a slope of ~132. Where the elliptic curve and the black line meet,
there is a cusp. Your elliptic curve doesn't have a cusp by itself, but the
joining of the elliptic curve and the continuation line does.
Perhaps I've misunderstood, and you don't care about a continuation leaving the
P1-P2 interval to the northeast. If so, then perhaps the elliptic curve meets
your needs.
If you look at the bottom graph, a repeat of an earlier one, you see the blue
cubic curve has a local minimum between P1 and P2. What Dave was saying is
this; If you want a curve that moves from P1 to P2, leaving P1 with a slope of
-0.1308565 and then approaching P2 with a slope of 132.5335679, you have two
choices of how to approach P2. You can approach with a northeasterly direction
of travel or with a southwesterly direction of travel. If you approach heading
to the northeast, you must be approaching from the left and below. Since the
y-coordinate of P2 is less than that of P1, you must first go below P2's y
coordinate in order to then approach while travelling northeasterly. This means
that there MUST be a local minimum in the cubic or any other curve that would
approach in the manner I described. But your elliptic curve has no minimum OF
THAT SORT. It does have a minimum where it meets P2, but it's the end of the
elliptic curve, not between P1 and P2 as is the local minimum of the blue cubic
curve.
Perhaps all of this has no relevance to what you want to do; perhaps what
happens outside of the interval P1 < x < P2 doesn't matter to you. But I
thought you wanted your curve to be a continuation in position and slope of the
data outside the interval.
I didn't mean (and I don't think Dave did either) that the cusp and minimum were
part of your elliptic curve. The cusp is due to the connection of the elliptic
curve with a path leading out of the interval to the northeast.
However, it's certainly possible to have a path leading away from P2 to the
southwest as a continuation of your elliptic curve after it reaches P2. It
would be an extension of the black line in the southwesterly direction. In that
case there will be no cusp, but I didn't think you would want to be following a
path that would lead your x-coordinate back into the interval between P1 and P2.
monir
>
> - Show quoted text -
Rodger;
Tremendous effort! Greatly appreciated and acknowledged!
1) >"I assumed that as you approach P2 from the left, you would want
to be moving up and to the right."
>"Imagine that you start at P1 and move on some curve toward P2. When you arrive at P2 I assume that you want to be moving on a slope of ~132, upward and to the right, NOT downward and to the left."
>"As I say, I have assumed that you want to be going northeast when you reach P2."
Why did you assume that ??
2) In an earlier reply, I summarized the problem and the conditions/
requirements of the sought function/curve joining P1(x1,y1) and
P2(x2,y2) as follows:
General:
x1 < x2 < 0
y1 > y2 > 0
sign of d1::(dy/dx) at P1 and d2::(dy/dx) at P2 could be any of the 4
+/- combinations.
Desired Solution:
a. satisfies the position and slope at the two endpoints;
b- numerically stable and does not require extra precision;
c. exact, smooth, continuous and well-defined;
d. interpolated y-value is bounded by y2 and y1 values; and
e. presents no extremum or inflection points within the range-of-
interest.
What happens to the interpolation curve leaving P2 is not and should
not be a concern.
3) >"I didn't mean (and I don't think Dave did either) that the cusp
and minimum were part of your elliptic curve. The cusp is due to the
connection of the elliptic curve with a path leading out of the
interval to the northeast."
I just wonder: Are you using a different definition of "cusp" ?? It
usually refers to a spike in the numerical values! In aeronautical
engineering, we say the trailing edge is "cusped" if the angle between
the upper and lower surfaces of the aerofoil is zero, a definition
which is analogous/consistent with the NA definition.
4) If you still see a "cusp due to the connection of the elliptic
curve with the path leading out ", please provide the corresponding x
value and the value of the "cusp".
5) I'm sorry but I don't see how the characteristics of the
"analytical" curve approaching P1 from the left and the curve leaving
P2 to the left (for the worked example) affect the elliptic solution
or produce a cusp.
Let us assume for the sake of discussion that both "analytical" curves
are straight-line segments with slope d1 and d2 respectively.
6) Should there be a better function/curve that satisfies a-d above, I
would be anxious to try it.
Any thoughts ?? Thanks again for your help and interest in the
problem. I find the discussion very useful!
Monir
The Phantom
<SNIP>
>Rodger;
>
>Tremendous effort! Greatly appreciated and acknowledged!
>
>1) >"I assumed that as you approach P2 from the left, you would want
>to be moving up and to the right."
>>"Imagine that you start at P1 and move on some curve toward P2. When you arrive at P2 I assume that you want to be moving on a slope of ~132, upward and to the right, NOT downward and to the left."
>>"As I say, I have assumed that you want to be going northeast when you reach P2."
>
>Why did you assume that ??
In an early post you said:
"I've a tiny "gap" in the analytical derivation between two points P1 and P2."
but then you say that you want to find a curve that has the proper derivative AT
the points P1 and P2. This makes it sound to me like the "gap" IS the entire
interval between P1 and P2. If the ends of the gap were not AT P1 and P2, you
would care about the derivative at points internal to P1 and P2 where the gap
begins and ends, not at exactly P1 and P2. And if the entire interval between
P1 and P2 IS, in fact, the "gap", then it would seem reasonable to assume that
there is something OUTSIDE the gap, otherwise it wouldn't be a "gap". The very
word "gap" implies something outside the "gap".
It would be EXTREMELY helpful if you would show, in a graph, the analytical
derivation you referred to, the "gap" in that derivation, and show P1 and P2 on
the plot.
>
>2) In an earlier reply, I summarized the problem and the conditions/
>requirements of the sought function/curve joining P1(x1,y1) and
>P2(x2,y2) as follows:
>General:
>x1 < x2 < 0
>y1 > y2 > 0
>sign of d1::(dy/dx) at P1 and d2::(dy/dx) at P2 could be any of the 4
>+/- combinations.
>Desired Solution:
>a. satisfies the position and slope at the two endpoints;
>b- numerically stable and does not require extra precision;
>c. exact, smooth, continuous and well-defined;
>d. interpolated y-value is bounded by y2 and y1 values; and
>e. presents no extremum or inflection points within the range-of-
>interest.
>What happens to the interpolation curve leaving P2 is not and should
>not be a concern.
"What happens to the interpolation curve leaving P2 is not and should not be a
concern."
This point has not been clear to me all along. I explain above why I assumed
that there was something "outside" the P1-P2 interval, for which you were trying
create a continuation inside the P1-P2 interval.
>
>3) >"I didn't mean (and I don't think Dave did either) that the cusp
>and minimum were part of your elliptic curve. The cusp is due to the
>connection of the elliptic curve with a path leading out of the
>interval to the northeast."
>
>I just wonder: Are you using a different definition of "cusp" ??
From my Dictionary of Mathematics:
"Cusp, A double point at which the two tangents to the curve are coincident. A
cusp of the first kind is a cusp in which there is a branch of the curve on EACH
side of the double tangent in the neighborhood of the point of tangency. A cusp
of the second kind is a cusp for which the two branches of the curve lie on the
SAME side of the tangent in the neighborhood of the point of tangency"
In this case the two tangents to the curve (the "curve" being the combination of
the elliptic curve, and the black continuation line) are the tangent to the
elliptic curve and the tangent to whatever continuation curve there may be. If
the continuation curve has the same slope at P2 as the elliptic curve, then they
form a cusp at P2 IF the continuation leaves P2 in a northeasterly direction.
Which kind of cusp depends on whether the continuation curve is concave up or
concave down. If the continuation curve is a straight line, then I guess it's a
cusp of the 1 and 1/2 kind. :-)
The dictionary gives an example: The semicubical parabola, y^2 = x^3, has a
cusp at the origin.
>It >usually refers to a spike in the numerical values!
I've been using it not to refer to a spike in numerical values, but in the
characteristic of a curve at a point, in accordance with the mathematical
definition given above.
Also, see:
http://en.wikipedia.org/wiki/Cusp_%28singularity%29
>In aeronautical
>engineering, we say the trailing edge is "cusped" if the angle between
>the upper and lower surfaces of the aerofoil is zero, a definition
>which is analogous/consistent with the NA definition.
Have a look at the graphs in:
http://img2.freeimagehosting.net/image.php?04b4672dfb.gif
The upper graph shows your elliptic curve in red, and a line with a slope of
132.5335679 meeting the elliptic curve at P2, P2 has an x value of -0.0248588.
I don't know what you mean by "value of the cusp". The elliptic curve and the
black line meet at an angle of zero, which would fit the aeronautical
engineering definition of "cusp".
The lower graph shows a meeting of the elliptic curve and the black line without
a cusp.
Are you able to see the graphs? Do you not agree that the upper graph shows a
cusp, and the lower graph does not? The definition from my Dictionary of
Mathematics applies to the top graph, but not the bottom graph.
>
>4) If you still see a "cusp due to the connection of the elliptic
>curve with the path leading out ", please provide the corresponding x
>value and the value of the "cusp".
See just above.
>
>5) I'm sorry but I don't see how the characteristics of the
>"analytical" curve approaching P1 from the left and the curve leaving
>P2 to the left (for the worked example) affect the elliptic solution
>or produce a cusp.
>Let us assume for the sake of discussion that both "analytical" curves
>are straight-line segments with slope d1 and d2 respectively.
In that case, the "analytical" curve can leave the elliptic curve at P2 in two
ways. It can leave going to the northeast, or it can leave going to the
southwest. Either can satisfy the slope requirement, but one case produces a
cusp (the top graph at the URL given above) and the other does not (bottom
graph).
>
>6) Should there be a better function/curve that satisfies a-d above, I
>would be anxious to try it.
Better is in the eye of the beholder. The cubic curve discussed earlier meets
P2 travelling northeasterly. Your elliptic curve meets P2 travelling
southwesterly. Does the difference not matter to you? If not, then the
elliptic curve may be "better" for your purposes.
monir
>
> - Show quoted text -- Hide quoted text -
>
> - Show quoted text -
Rodger;
Clearly there has been some misunderstanding, and it was probably my
fault as I'll explain below.
1) First, I've been able to review your latest graphs in:
http://img2.freeimagehosting.net/image.php?04b4672dfb.gif
Here's a possible cause of the confusion!
2)>"The upper graph shows your elliptic curve in red, and a line with
a slope of 132.5335679 meeting the elliptic curve at P2."
>"The elliptic curve and the black line meet at an angle of zero, which would fit the aeronautical engineering definition of cusp". >"(the "curve" being the combination of the elliptic curve, and the black continuation line)"
>"Do you not agree that the upper graph shows a cusp, and the lower graph does not?
I conditionally agree! In the upper graph, the black line is the
"tangent line" at P2 and it is NOT a "continuation line". It is not
an analytical line nor part of the solution yet it appears to be the
cause of the confusion!
You can extend it to infinity on both directions if you choose to,
ignore it completely, or remove it all together. It wouldn't matter!
I decided earlier on to display the "slope line" at P2 just to show
that the elliptic curve satisfies the slope condition and smoothly
approaches the endpoint P2 with no cusps, not realizing that the
display would cause such confusion!!!!!!
3)>"I don't know what you mean by "value of the cusp".
I was under the impression all along that you were using cusp and
spike interchangeably! Your example of y = x^(2/3) is a good one for
zero cusp value at origin. If I'm not mistaken, some curves such as
the Tractrix family of curves have finite cusp values (I think!). In
any event, I suppose it's a mute point following the above
clarifications!
Thanks for the exact and clear definition from the Dictionary of
Mathematics.
4)>"In that case, the "analytical" curve can leave the elliptic curve
at P2 in two ways. It can leave going to the northeast, or it can
leave going to the southwest. Either can satisfy the slope
requirement, but one case produces a
cusp (the top graph at the URL given above) and the other does not
(bottom graph)."
The bottom graph represents the ONLY correct way for the worked
example:
- it satisfies all the conditions/requirements (items 2.a to 2.e of my
previous post);
- it doesn't display the apparently ever-confusing slope line at P2;
and
- it smoothly and nicely joins the elliptic curve with the (assumed)
black analytical straight-line segment at P2.
Does the above clear up the confusion ??
Regards.
Monir
The Phantom
<SNIP>
>
Somewhat. But the black "continuation" line in the bottom graph "re-enters" the
P1-P2 interval. It goes "backwards" in the x direction. Is that ok with you?
When you were considering the cubic curve, it meets P2 going in the northeast
direction, and if you were to combine it with the black line in the top graph,
there would be no cusp, because the direction of travel would just be
continuously northeast. At that time you didn't express dissatisfaction with
that aspect of the cubic.
I guess you don't care about "continuation" lines, and if the elliptic curve
meets your needs, so be it.
I think your use of the word "gap" confused me. I think now that you just
wanted to find a curve to connect P1 and P2 within some constraints, and by the
word "gap" you just meant that you hadn't found such a function. But I was led
to believe that there was some curve or function to the left of P1 and the right
of P2 that you wanted to continue in the P1-P2 interval. Apparently this was a
wrong impression.
I'm curious. Is this about designing airfoils?
>
>Regards.
>Monir
monir
>
> - Show quoted text -- Hide quoted text -
>
> - Show quoted text -
Rodger;
Good! It seems that the misunderstanding has been cleared and everyone
is happy now!
>I'm curious. Is this about designing airfoils?
It's primarily related to the study of fluid flow dynamics associated
with lifting surfaces and a continuation of my personal research/
teaching activities.
Thanks again for your help, patience and interest.
Regards.
Monir