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.