Sunday, August 30, 2020

Numworks: Quadratic Equation and Newton's Method - Do Guesses Matter?

Numworks: Quadratic Equation and Newton's Method - Do Guesses Matter?  

Introduction

Let A, B, C, be positive, non-zero coefficients of the quadratic equation:

Ax^2 + Bx - C = 0 with A > 0, B > 0, C > 0

First, start with an initial guess x0 = 0.  Count how many iterations it takes for the method to reach a suitable answer.   For this experiment, we will use a tolerance of 8 decimal places (tol = 1E-8), setting a maximum of 200 iterations.   

Consider a "reduced" equation to obtain guesses for the original equation:  

Ax^2 - C = 0

Hence, our initial guesses are x0 = √(C/A) and x0 = -√(C/A)

Using Newton's Method to determine the roots with a guess of x0:

x_n+1 = x_n - [Ax^2 + Bx - C]/[2Ax + B]

So which of the following guesses gets us to a solution faster?

x0 = 0

x0 = √(C/A) 

x0 = -√(C/A)


The Python script below tests the question.  

Python Test Script:  guesttest.py


from math import *

def testf(a,b,c,w0):

 n=1

 w1=1

 while abs(w1)>1e-8

  w1=(a*w0**2+b*w0-c)/(2*a*w0+b)

  w0=w0-w1

  n=n+1

  if n==201:

   print('Exceeded iterations')

   break

 print('Root: '+str(w0)

 print('Iterations: '+str(n))

 return


# EWS 2020-08-10

print("y=ax**2+bx-c=0")

print("a>0,b>0,c>0")

a=float(input('a? '))

b=float(input('b? '))

c=float(input('c? '))

# initial conditions

x1=0

x2=sqrt(c/a)

x3=-sqrt(c/a)

# results

print('x1: '+str(x1))

testf(a,b,c,x1)

print('----')

print('x2: '+str(x2))

testf(a,b,c,x2)

print('----')

print('x3: '+str(x3))

testf(a,b,c,x3)

print('----')


Results

Test Results - Iterations Used to get to a Solution

Ax^2 + Bx - C = 0 with A > 0, B > 0, C > 0


A B C x0 = 0 x0 = √(C/A) x0 =-√(C/A)
1 1 4 6 6 8
1 4 4 6 6 error reached
1 5 4 6 6 9
1 5 6 6 6 12
2 5 6 6 6 8
3 5 6 7 6 7
8 13 81 8 6 6
16 13 81 8 5 6
32 13 81 9 5 5
64 13 81 9 5 5
1 40 10 5 6 6
1 40 25 5 6 6
1 40 40 5 6 7
1 40 55 6 6 8
10 18 42 7 6 7
10 18 52 7 6 6
10 18 62 7 6 6
10 18 72 7 6 6

From the data presented, I think it is inconclusive of whether which guess generally makes for getting to the solution faster.   What do you find?  

Eddie 

All original content copyright, © 2011-2020.  Edward Shore.   Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited.  This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author. 






Saturday, August 29, 2020

Radio Shack EC-4019: Catalan and Borel Triangles

 Radio Shack EC-4019:   Catalan and Borel Triangles 

Introduction:   Catalan and Borel Triangles 

Catalan numbers counts a number of lattice paths, straight lines with segment of length 1, in the Cartesian plane from the origin (0,0) to the point (n,k).  There path does not go above the line y = x.  

Entries in Borel's Triangle depend on the entries of Catalan's Triangle.  There are various interpretations of how counting relate to Borel's triangle, including counting the number of Dyck path of semi-length n+1, counting path in leaf-market binary trees, and set of binary trees with n+1 vertices with k marked vertices with no markings on the tree's right spine.  Please see the source below for more details.  


Formula Derivation

You can find any entry in Catalan's triangle by the formula:

C = (n - k + 1) / (n + 1) * COMB(n + k, n)

where:


1.  k ≤ n, I do believe that k and n need to be integers, and

2.  COMB is the nCr, combination function, where for any x and y:

COMB(x, y) = x! / ( y! * (x - y)! )

We can simplify the formula for Catalan's triangle by:

C 

= (n - k + 1) / (n + 1) * COMB(n + k, n)

= (n - k + 1) / (n + 1) * (n + k)! / (n! * (n + k - n)!)

= (n - k + 1) / (n + 1) * (n + k)! / (n! * k!)

= (n - k + 1) / (n + 1) * (n + k)! /n! * 1/k!


Since (n + 1) * n! = (n + 1) * n * (n - 1) * (n - 2) * ... = (n + 1)!

= (n - k + 1) * (n + k)! / ((n + 1)! * (n - k)!)


Observe that

C / (n-k+1) = (n + k)! / ((n + 1)! * (n - k)!)


Similarly, you can find any entry in Borel's triangle by the formula:


B 

= 1 / (n + 1) * COMB(2n + 2, n - k) * COMB(n + k, n)

= 1/ (n + 1) * (2n + 2)! / ((n - k)! * (n + 2 + k)!) * (n + k)!/(n! * k!)

= (2n + 2)! / ((n - k)! * (n + 2 + k)!) * (n + k)!/((n+1)! * k!)

= (2n + 2)! / ((n - k)! * (n + k + 2)!) * Ca


Radio Shack EC-4019 Programs:  Catalan and Borel Triangle Numbers


To run the program:

1.  Store k in memory register 1.  Keystrokes:  k [ Kin ]  [ 1 ]

2.  Store n in memory register 2.  Keystrokes:  n [ Kin ] [ 2 ]

3.  To find the Catalan Triangle, run program I.   

4.  To find the Borel triangle number, run program I, then immediately run program II.  Program II depends on the result of program I.  


Make sure that k ≤ n, as the programs assume that your inputs are valid.


Program I:  Catalan Triangle Number

(25 steps)

(small x:  multiply key,  slash: divide key)


[(---

Kout 2

-

Kout 1

+


1 

---)]

x 

[(---

Kout 2


+

Kout 1

---)]

x!

/


[(---

Kout 2

+

1

---)]


x!

/

Kout 1

x!

=


Program II: Borel Triangle Number

(28 steps)

(small x:  multiply key,  slash: divide key)


x

[(---

2

x 

Kout 2


+

2

---)]

x!

/


[(

Kout 2

-

Kout 1

+


1

)]

x!

/

[(


Kout 1

+

Kout 2

+

2


---)]

x!

=


Examples

k = 2,  n = 6,  Catalan (I):  20,  Borel (II): 4004

k = 3, n = 5,  Catalan (I): 28, Borel (II): 616

k = 4, n = 4,  Catalan (I): 14, Borel (II): 14


Source:

Cai, Yue and Yan, Catherine.  "Coutning with Borel's Triangle"  Elsevier B.V. Discrete Mathematics.   November 15, 2018

https://arxiv.org/pdf/1804.01597.pdf

Also:  https://doi.org/10.1016/j.disc.2018.10.031

Retrieved August 9, 2020

Eddie

All original content copyright, © 2011-2020.  Edward Shore.   Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited.  This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author. 


Sunday, August 23, 2020

Using Grand Total and Independent Memory for Complex Invoicing

Using Grand Total and Independent Memory for Complex Invoicing 

Introduction

Portable desktop (usually solar powered) calculators usually have at least two memory registers:  independent memory (M) and grand total (GT).  

Adding calculations to the independent memory by pressing the memory plus key    [ M+ ].  

Adding calculations to the grand total by pressing the equals key [ = ]. 

Invoice with Taxable and Nontaxable Items

We have an invoice that has items that are subject to sales tax, and items and services that are exempt from sales tax. 

Take advantage of both the independent and grand total memories to calculate the total invoice.

Two approaches:

[ ON/C ] to clear the grand total memory
[ MC ] (or press [ MRC ] until the independent memory is cleared) 

Method 1: 
taxable items [ M+ ]
nontaxable items [ = ]

After finishing entering all the items, add sales tax:
[ MR ] [ × ] sales tax rate [ % ] [ M+ ]

Calculate the total invoice:
[ MR ] [ = ] [ GRAND TOTAL ]

Method 2:
taxable items [ = ]
nontaxable items [ M+ ]

After finishing entering all the items, add sales tax:
[ GRAND TOTAL ] [ × ] sales tax rate [ % ] [ = ]

Calculate the total invoice:
[ MR ] [ = ] [ GRAND TOTAL ]

The example will use Method 1.

An Example of an Invoice

Taxable Items:
$ 19.99
$ 39.96
$ 14.97

Nontaxable Items:  
$ 109.00
$  15.00

Sales Tax Rate:  9% 

As long as you keep your designation of what memory registers are used, you can enter each item in any order.  So both methods should work.

[ AC/ON ] [ MC ]
109 [ = ]
15 [ = ]
19.99 [ M+ ] 
39.96 [ M+ ]
14.97 [ M+ ]
[ MR ] [ × ] 9 [ % ] [ M+ ] 
[ MR ] [ = ] 
[ GRAND TOTAL ]

Result:  205.6628

[ AC/ON ] [ MC ]
109 [ = ]
19.99 [ M+ ] 
15 [ = ]
39.96 [ M+ ]
14.97 [ M+ ]
[ MR ] [ × ] 9 [ % ] [ M+ ] 
[ MR ] [ = ] 
[ GRAND TOTAL ]

Result:  205.6628

Invoice:  $205.66

The calculators I used for this example are the Casio JF-100BM and Casio MC-12M Shop Calculator.  

Eddie

All original content copyright, © 2011-2020.  Edward Shore.   Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited.  This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author. 

Saturday, August 22, 2020

HP Prime: Psuedo Ordering of Complex Numbers

HP Prime: Psuedo Ordering of Complex Numbers 


Can Complex Numbers be Ordered?

Not really, not if you consider what is required to be an ordered field.  The complex numbers can not be an ordered field, one counter example is illustrating the nonzero complex numbers i and 1, when squared and added together, you get zero.

i^2 + 1^2 = -1 + 1 = 0

If you order complex numbers by their modulus (absolute values) only, you run into the problem of having many complex numbers becoming equivalent.  Example:

2 + 3i  and 3 + 2i,  each with radius √13

One way to "psuedo" order complex numbers is to use the numbers looking at both their modulus and angle.   

z1  "<"  z2 if and only if one of the two conditions are met:

| z1 | < | z2 | or

| z1 | = | z2 | and atan(imag(z1)/real(z1)) < atan(imag(z2)/real(z2))

I believe we do not have about the unit of angle  (degrees vs. radians vs. grads) for this psuedo-ordering.

Please see the sources listed at the end of this article for more details. 


HP Prime Program:  CLESS

EXPORT CLESS(z1,z2)

BEGIN

// 2020-08-06 EWS

// psuedo ordering of 

// complex numbers

LOCAL a1,a2,t1,t2;

a1:=ABS(z1);

a2:=ABS(z2);

t1:=ATAN(IM(z1)/RE(z1));

t2:=ATAN(IM(z2)/RE(z2));

IF (a1<a2) OR ((a1==a2) AND (t1<t2)) THEN

RETURN "TRUE";

ELSE 

RETURN "FALSE";

END;

END;


Example (angle shown in radians, rounded to five decimal places):


These complex numbers are listed in order:

1 + i    (1.41421 ∠ 0.78540)

2 + i    (2.23067 ∠ 0.46365)

1+ 2i  (2.23607 ∠ 1.10715)

2 + 2i  (2.82843 ∠ 0.78540)

3 + i  (3.16228 ∠ 0.32175)

1 + 3i (3.16228 ∠ 1.24905)

3 + 2i (3.60555 ∠ 0.58800)

2 + 3i (3.60555 ∠ 0.98279)

4 + i (4.12311 ∠ 0.24498)


Sources:

Weimer, Richard C.   "Can the Complex Numbers Be Ordered?"  The Two-Year College Mathematics Journal, Dec. 1976, Vol. 7, No. 4.   Taylor & Francis, Ltd on behalf of the Mathematical Association of America.  https://www.jstor.org/stable/3027050?seq=1


Yada, Dharmendra Kumar.  "A New Approach to Ordering Complex Numbers"  International Journal of Mathematical Sciences and Engineering Applications (IJMSEA), Vol 2.  No. III (2008), pp. 221-223.   Article downloaded from ResearchGate:  https://www.researchgate.net/publication/267465398_A_new_approach_to_ordering_complex_numbers

Eddie

All original content copyright, © 2011-2020.  Edward Shore.   Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited.  This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author. 


Sunday, August 16, 2020

Review: Casio MC-12M Shop Calculator

Review:  Casio  MC-12M Shop Calculator



Just The Facts:

Model:  MC-12M Shop Calculator 
Company: Casio
Type:  Desktop 
Display:  12 digits
Battery:  Solar and Battery (AA battery backup)
Logic: Chain
Memory:  Two:  M and Grand Total 
Double Zero Key


I purchased this calculator several months ago at a WalMart for about $11.50. 

Features:

Cost/Sell/Margin Calculation
Change Function
Square Root

The display is set an angle so it is tilted towards the user.  I find the keyboard a pleasure to use.  My hands may be a little big, but with the right size of hands you can speed type on the MC-12M.   The keys are fast and responsive.  

Memory

There are two memories within the MC-12M:  the independent memory M and the Grand Total key.  

The independent memory works with the usual memory plus (M+) and memory minus (M-) keys.   The recall (MR) and clear memory (MC) keys are separate keys on the MC-12M.  

The grand total memory is activated every time you press the equals [ = ] key.  On the MC-12M, the grand total memory is automatic.   The grand total memory is cleared when either [AC/ON] is pressed or the calculator is turned off.  

The Change Feature

The big selling feature of the MC-12M Shop Calculator is the [CHANGE] button.   You enter an amount and press [CHANGE], like you were handed cash in a transaction.   When it is pressed, it compares the amount entered against the amount in the grand total memory registered.  The independent memory is not used in the change feature.  

Example:

You have a $100 bill (C-note) and purchase the following:
1 SD card for $19.95
1 desktop calculator for $11.49  ;)
3 reams of paper for $5.00 each (on sale)
15 pens at $0.49 each
2 USB drives, 32GB for $10.60 each
Sales Tax: 9%

Keystrokes:

[AC/ON]
19.95 [ = ]
11.49 [ = ]
3 [ × ] 5 [ = ]
15 [ × ] .49 [ = ]
2 [ × ] 10.60 [ = ]
[ GRAND TOTAL ]     (Display: 74.99)

[ × ] 9 [ % ] [ = ]

[ GRAND TOTAL ]     (Display/Total Invoice: 81.7391   [$81.74]  )

100 [ CHANGE ]   (Change: 18.2609  [$18.26] )

Percent Change

You can calculate percent change with the use of the minus and percent keys.

%CHG Keystrokes:   new [ - ] old [ % ]

Calculate the percent change from 14.75 to 19.85?

19.85 [ - ] 14.75 [ % ] 

Percent Change:  34.5762711864   (34.5762711864%)

Wishes and Verdict

There are two things the MC-12M Shop Calculator have:

Fix Decimal Selector:  This would be great for the shop because most currency calculations work with 2 decimal places.   My selection of choices would be 0, 1, 2, 3, 4, 6.

I would also like to see the Add Mode added. 

Since this is a shop calculator, I can imagine people using this calculator to calculator areas and volumes, so a Ï€ key would be nice.  

Ï€ ≈ 3.1415926535

This calculator is nice to operate and hold and yes, I like the square root (which should be present on more desktop calculators).   Worth a look.

Eddie

All original content copyright, © 2011-2020.  Edward Shore.   Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited.  This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author. 

Saturday, August 15, 2020

TI 84 Plus CE: Cauchy Fit Regression

TI 84 Plus CE: Cauchy Fit Regression

Note, in this blog capital and lower case letter variables are not assumed to be the same,  hence A ≠ a, B ≠ b, C ≠ c.

Derivation

General Equation:

Y = 1 / (A * (x + B)^2 + C)

We can set the equation in a general quadratic equation format by:

Y = 1 / (A * (x + B)^2 + C)
1/Y = A * (x + B)^2 + C
1/Y = A * (x^2 + 2*B*x + B^2) + C
1/Y = A*x^2 + 2*A*B*x + A*B^2 + C
1/Y = A*x^2 + 2*A*B*x + [A*B^2 + C]

Let:
y = 1/Y
a = A
b = 2*A*B
c = A*B^2 + C

Then:
y = a*x^2 + b*x + c

This form can be used when using a graphing or scientific calculator to calculate quadratic regression. 

We can solve for A, B, C:

A = a

b = 2*A*B
B = b / (2*A)

c = A*B^2 + C
C = c - A*B^2

To find a Cauchy regression fit:

1.  Enter the x data.
2.  Enter the y data but take the reciprocal of each y data point. 
3.  Run the Quadratic Regression:  y = a*x^2 + b*x + c
4.  Solve for the Cauchy regression to get the coefficients:

A = a
B = b / (2*A)
C = c - A*B^2

TI-84 Plus CE Program:  CAUCHFIT
Size:  200 Bytes

"2020-07-31 EWS"
Disp "CAUCHY FIT"
Input "X DATA? ", L1
Input "Y DATA? ", L3
1 / L3 → L2
QuadReg L1, L2
a → A
b / (2 A) → B
c - A B^2 → C
ClrHome
Output(3,1,"Y=1/(A(X+B)^2+C)"
Output(4,1,"A=")
Output(4,3,A)
Output(5,1,"B=")
Output(5,3,B)
Output(6,1,"C=")
Output(6,3,C)
Pause

Notes:

1.  The variables a, b, and c are from the VARS, Statistics, EQ menu.

2.  x data is stored in L1, 1/y data is stored in L2, and y data is stored in L3.  To call lists L1, L2, and L3, press [ 2nd ] [ 1 ], [ 2nd ] [ 2 ], and [ 2nd ] [ 3 ], respectively. 

3.  The use of Output should allow this code to be programmed on earlier TI-83 Plus and TI-84 Plus calculators. 

Example

X = { 1, 1.5, 2, 2.5 , 3}
Y = { 0.0774, 0.0577, 0.0442, 0.0347, 0.0278}

Results:

A = 1.823972632
B = 1.157397925
C = 4.437539736

y ≈ 1 / ( 1.823972632 * (x + 1.157397925)^2 + 4.437539736 )

Source:

Kolb, William M.  Curve Fitting For Programmable Calculators IMTEC.  Bowie, MD 20716.  1982.  ISBN-10:  0-943494-00-01


Stay safe and sane everyone.  Blessings,

Eddie

All original content copyright, © 2011-2020.  Edward Shore.   Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited.  This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author.

Sunday, August 9, 2020

HP 12C: Last X: Sums and Products

HP 12C:   Last X:   Sums and Products

Introduction

It is possible to use the LST X feature on the HP 12C calculators assist us in quickly calculate sums or products  of a terms, especially when the list of terms are in a series. 

For today's blog, we are starting with a base amount, a, and then adding 1 to each term.   You can use a similar algorithm for a sequence where each term is doubled, tripled, 1 is subtracted from the previous term, and so on.  The key is to complete the adjustment, use storage arithmetic, and then use LST X.

Accessing LST X:

HP 12C (Classic):  [ g ] [ ENTER ]

HP 12C Platinum:  [ g ] [ + ]

These algorithms can be used in program. 


Sum:  a + (a+1) + (a+2)  + (a+3) + ....

Let n represent the storage register to be used.  On the HP 12C, only storage registers R0 through R4 (classic HP 12C).

Algorithm:

a 
STO n

Loop:  
LST X 
1     [ or 2, x  to double each term,  1, - to subtract 1, etc.]
+  
STO+ n    [stores what is in the X display to the Last X register]

Finish:  
RCL n


Example:  7 + 8 + 9 +10 = 34.   Use register 0 to store the sum.

7
STO 0

LST X
1
+
STO+ 0

LST X
1
+
STO+ 0

LST X
1
+
STO+ 0

RCL 0


Sum:  1/a + 1/(a+1) + 1/(a+2) + 1/(a+3) + ....

Algorithm:

a 
1/x
STO n

Loop:  
LST X 
1     
+  
1/x
STO+ n    

Finish:  
RCL n

Example:  1/7 + 1/8 + 1/9 + 1/10 ≈ 0.47897

7
1/x
STO 0

LST X
1
+
1/x
STO+ 0

LST X
1
+
1/x
STO+ 0

LST X
1
+
1/x
STO+ 0

RCL 0

Sum:  √a + √(a+1) + √(a+2) + √(a+3) + ...

Algorithm:

a 
√
STO n

Loop:  
LST X 
1     
+  
√
STO+ n    

Finish:  
RCL n

Try the algorithm on this example:  √7 + √8 + √9 + √10 ≈ 11.63646

Sum:  a^2 + (a+1)^2 + (a+2)^2 + (a+3)^2 + ...

Hint:  Use ENTER, x instead of 2, y^x

Algorithm:

a 
ENTER 
*
STO n

Loop:  
LST X 
1     
+  
ENTER
*
STO+ n    

Finish:  
RCL n

Try the algorithm on this example:  7^2 + 8^2 + 9^2 + 10^2 = 294

Let's move on to products.

Product:  a * (a+1) * (a+2) * (a+3) * ....

a 
STO n

Loop:  
LST X 
1     
+  
STOx n    

Finish:  
RCL n


Example:  7 * 8 * 9 * 10 = 5040.   Use register 0 to store the product.

7
STO 0

LST X
1
+
STOx 0

LST X
1
+
STOx 0

LST X
1
+
STOx 0

RCL 0

Product:  1/a * 1/(a+1) * 1/(a+2) * 1/(a+3) * ....

a 
1/x
STO n

Loop:  
LST X 
1     
+  
1/x
STOx n    

Finish:  
RCL n


Try the algorithm on this example 1/7 * 1/8 * 1/9 * 1/10 ≈ 0.0020

Product:  √a * √(a+1) * √(a+2) * √(a+3) * .......

a 
√
STO n

Loop:  
LST X 
1     
+  
√
STOx n    

Finish:  
RCL n


Try the algorithm on this example √7 * √8 * √9 * √10 ≈ 70.99296

Hopefully this will help your calculations on your RPN calculators faster and more effective,

Eddie

All original content copyright, © 2011-2020.  Edward Shore.   Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited.  This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author.

Saturday, August 8, 2020

TI-84 Plus CE: Random Numbers and Linear Regression

Random Numbers and Linear Regression


The Test Case y = rand + x

The goal is to see how good we can fit a line to the points generated by the equation:

y = rand + x

where rand is a random generated number between 0 + 1.  Plotting such equation will show that y = rand + x will be bounded in between y = x and y = x + 1 (see screen shot below):




In each of the three tests cases presented, let x be a range of integers between 0 and 14. 

X = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14}

In the cases presented today, all random amounts are rounded to three decimal places.

Case 1

Y = {0.735, 1.931, 2.007, 3.029, 4.499, 5.578, 6.794, 7.235, 8.847, 9.485, 10.326, 11.873, 12.287, 13.795, 14.101}

Intercept ≈ 0.530
Slope ≈ 0.996
Correlation ≈ 0.997

Case 2

Y = {0.900, 1.792, 2.462, 3.639, 4.083, 5.869, 6.088, 7.225, 8.458, 9.394, 10.401, 11.329, 12.283, 13.441, 14.478}

Intercept ≈ 0.627
Slope ≈ 0.976
Correlation ≈ 0.999

Case 3

Y = {0.066, 1.687, 2.654, 3.166, 4.879, 5.098, 6.891, 7.235, 8.504, 9.615, 10.217, 11.040, 12.823, 13.725, 14.397}

Intercept ≈ 0.428
Slope ≈ 1.006
Correlation ≈ 0.998

Testing 100 Curve Fits

The program RANDLR will test 100 test cases of 50 data points,  X is a list generated of the integers from -24 to 25, while Y is calculated by X + rand, rounded to three decimal places.  The intercept, slope, and correlation shown are the arithmetic average of 100 test runs.   I used various random seeds to here are the results I used from five test runs:

Test 1:
Intercept ≈ 0.4977
Slope ≈ 1.0007
Correlation ≈ 0.9998

Test 2:
Intercept ≈ 0.4976
Slope ≈ 1.0007
Correlation ≈ 0.9998

Test 3:
Intercept ≈ 0.5012
Slope ≈ 1.0002
Correlation ≈ 0.9998

Test 4:
Intercept ≈ 0.5046
Slope ≈ 0.9997
Correlation ≈ 0.9998

Test 5:
Intercept ≈ 0.5012
Slope ≈ 1.0002
Correlation ≈ 0.9998

I can roughly estimate that a line to estimate y = x + rand is y = x + 0.5.

TI-84 Plus CE Program:  RANDLR
Size:  322 bytes

"2020-07-26 EWS"
ClrHome
Disp "TEST:  Y=AX+B","AGAINST Y=X+rand","PLEASE WAIT."
Fix 4
seq(X,X,-24,25,1) → L1
0 → A
0 → B
0 → R
For(I,1,100)
Output(5,10,"      ")    \\ six spaces
Output(5,1,"PROGRESS:")
Output(5,12,I)
seq(round(X+rand,3),X,-24,25,1) → L2
LinReg(ax+b) L1, L2
A+a → A       
B+b → B
R+r → R
End
.01 A → A
.01 B → B
.01 R → R
ClrHome
Output(3,1,"AVERAGE")
Output(4,1,"ITC:")
Output(4,6,B)
Output(5,1,"SLP:")
Output(5,6,A)
Output(6,1,"COR:")
Output(6,6,R)

Note: 

The lines

Output(5,10,"      ")    \\ six spaces
Output(5,1,"PROGRESS:")
Output(5,12,I)

generate an "in progress" counter so that the user how much longer the program will take.

L1 represents list 1 by pressing [ 2nd ] [ 1 ]. Likewise, L2 stands for list 2 and you can get the L2 character by pressing [ 2nd ] [ 2 ]. 

The statistics variables a, b, r are found in the VARS, Statistics, EQ menu.

The code above should work on any of the TI-83 Plus/TI-84 Plus family, but I only used it on the current TI-84 Plus CE. 

Until next time,

Eddie

All original content copyright, © 2011-2020.  Edward Shore.   Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited.  This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author.

Sunday, August 2, 2020

HP Prime and Casio fx-CG 50: Binary Conversions

HP Prime and Casio fx-CG 50:   Binary Conversions 

x = 2 * (1 + n)^p

The program DECOM2 takes a real number, x, and decomposes it to the form of:

x = 2 * (1 + n)^p

HP Prime Program:  DECOMP2

EXPORT DECOMP2(x)
BEGIN
// 2020-07-11 EWS
// x → 2^p*(1+n)
LOCAL s,p,n,str;
s:=ABS(x);
p:=FLOOR(LOG(s)/LOG(2));
n:=s/2^p-1;
str:=STRING(SIGN(X))+"*2^"+STRING(p)+
"*(1+"+STRING(n)+")";
RETURN str;
END;

Casio fx-CG 50 Program:  DECOMP2
188 bytes

"2020-07-11 EWS"
"X"? → X
Abs X → S
Int (log S ÷ log 2) → P
S < 1 ⇒ P - 1 → P 
S ÷ 2^P - 1 → N
X ÷ S → T
ClrText
Locate 1, 3, X
Locate 18, 3, "="
Locate 1, 5, T
Locate 4, 5, "×2^"
Locate 9, 5, P
Locate 18, 5, "×"
Locate 1, 6, "(1+"
Locate 5, 6, N
Locate 18, 6, ")"

Examples:

-720 = -1 * 2^9 * (1 + 0.40625)

0.7868 = 1 * 2^-1 * (1 + 0.5736)

Converting Numbers to Binary

The program RBIN converts any real number in base 10 to base 2, including non-integers and negative numbers.  The result is displayed as a string.   The HP Prime program has one string for output while the Casio fx-CG 50 uses two:  Str1 for the integer part and Str2 for the fractional part.

RBIN will need two arguments:  the number to be converted and the precision level.

For the precision level:

1 rounds the number to the nearest 1/2 before conversion.

2 rounds the number to the nearest 1/4 before conversion.

3 rounds the number to the nearest 1/8 before conversion.

4 rounds the number to the nearest 1/16 before conversion.

and so on.

HP Prime Program:  RBIN

EXPORT RBIN(x,p)
BEGIN
// 2020-07-11 EWS
// Decimal to Binary
// decimal, precision
// for real x, rounded 1/2^p

LOCAL b,i,f,n,s,w,k,str;
s:=ABS(x);
b:=IP(s)+ROUND(FP(s)*2^p,0)/2^p;
str:="";
i:=IP(b);
f:=FP(b);

IF i≠0 THEN
n:=IP(LOG(i)/LOG(2));
FOR k FROM n DOWNTO 0 STEP 1 DO
w:=IP(i/2^k);
str:=str+STRING(w);
i:=i-w*2^k; 
END;
END;

IF f≠0 THEN
str:=str+".";
FOR k FROM 0 TO p-1 DO
w:=IP(2*f);
str:=str+STRING(w);
f:=FP(2*f);
END;
END;

IF x<0 font="" then="">
str:="-"+str;
END;

RETURN str;
END;

Casio fx-CG 50:  RBIN
372 bytes

"2020-07-11 EWS"
"DEC>BIN: NUM"? → X
"PRECISION"? → P

Abs X → S
Int S + RndFix(Frac S × 2^P, 0) ÷ 2^P → B
" " → Str 1
Int B → I
Frac B → F

If I ≠ 0
Then
Int (log I ÷ log 2) → N
For N → K To 0 Step -1
Int (I ÷ 2^K) → W
W = 0 ⇒ Str 1 + "0" → Str 1
W = 1 ⇒ Str 1 + "1" → Str 1
I - W × 2^K → I
Next
IfEnd

" " → Str 2
If F ≠ 0
Then
"." → Str 2
For 0 → K To P-1
Int (2F) → W
W = 0 ⇒ Str 2 + "0" → Str 2
W = 1 ⇒ Str 2 + "1" → Str 2
Frac (2F) → F
Next
IfEnd

X < 0 ⇒ "-" + Str 1 → Str 1

ClrText
Locate 1, 3, X
Locate 20, 3, "="
Locate 1, 5, Str 1
Locate 1, 6, Str 2

Examples:

RBIN(-0.985,3) returns "-1"
RBIN(-0.985,6) returns "-.111111"
RBIN(-0.985,9) returns "-.111111000"

RBIN(860.63,3) returns "1101011100.101"
RBIN(860.63,6) returns "1101011100.101000"
RBIN(860.63,9) returns "1101011100.101000011"

Take care,

Eddie

All original content copyright, © 2011-2020.  Edward Shore.   Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited.  This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author.

Saturday, August 1, 2020

Casio fx-9750GIII and fx-CG50: System of Two Differential Equations, Runge Kutta 4th Order

Casio fx-9750GIII and fx-CG50: System of Two Differential Equations, Runge Kutta 4th Order

Introduction




The program TWORK4 uses the Runge Kutta 4th Order method to solve the following system of differential equations:

dx/dt = f(t, x, y)
dy/dt = g(t, x, y)

with initial conditions x0 = x(t0) and y0 = y(t0)

The next step is calculated with step h from:

x1 ≈ x0 + (k1 + 2 * k2 + 2 * k3 + k4) / 6
y1 ≈ y0 + (l1 + 2 * l2 + 2 * l3 + l4) / 6

k1 = h * f(t0, x0, y0)
l1 = h * g(t0, x0, y0)

k2 = h * f(t0 + h / 2, x0 + k1 / 2, y0 + l1 / 2)
l2 = h * g(t0 + h / 2, x0 + k1 / 2, y0 + l1 / 2)

k3 = h * f(t0 + h / 2, x0 + k2 / 2, y0 + l2 / 2)
l3 = h * g(t0 + h / 2, x0 + k2 / 2, y0 + l2 / 2)

k4 = h * f(t0 + h, x0 + k3, y0 + l3)
l4 = h * g(t0 + h, x0 + k3, y0 + l3)

For the next step set t0 = t0 + h, x0 = x1, and y0 = y1. 

Inputs:

DX/DT:  Enter dx/dt as a function of T, X, and Y.   T is the independent variable.

DY/DT:  Enter dy/dt as a function of T, X, and Y.   T is the independent variable.

T0, X0, Y0:  Enter the initial conditions

STEP:  Enter the step size

ITERATIONS: Enter the number of iterations desired.  This allows you to calculate a far point in one leap.  Example, if your initial condition is to = 0 and you want to find the point when t = 1 using the step size of 0.1, enter 0.1 for STEP and 10 for ITERATIONS.

Casio fx-9750GIII and fx-CG 50 Program: TWORK4

Notes: 

*  Although the code for the two calculators are the same, the programs will need to be programmed on each calculator separately.

*  The slash character (/) is accessed from the CHAR submenu.  The submenu shows up at the top-level menu which would read:
TOP/BOTTOM/SEARCH/MENU/A ←→ a/CHAR

* fn1 and fn2 are stored as function memories. 

"2020-06-18 EWS"
Rad
"DX/DT="? → fn1
"DY/DT="? → fn2
"T0"? → U
"X0"? → A
"Y0"? → B
"STEP?" → H
4 → Dim List 25
4 → Dim List 26
Lbl 0
"ITERATIONS"? → N
For 1 → I To N
A → X: B → Y: U → T
H * fn1 → List 25[1]
H * fn2 → List 26[1]
A + List 25[1] ÷ 2 → X
B + List 26[1] ÷ 2 → Y
U + H ÷ 2 → T
H * fn1 → List 25[2]
H * fn2 → List 26[2]
A + List 25[2] ÷ 2 → X
B + List 26[2] ÷ 2 → Y
H * fn1 → List 25[3]
H * fn2 → List 26[3]
A + List 25[3] → X
B + List 26[3] → Y
U + H → T
H * fn1 → List 25[4]
H * fn2 → List 26[4]
A + (List 25[2] + List 25[3] + Sum List 25) ÷ 6 → A
B + (List 26[2] + List 26[3] + Sum List 26) ÷ 6 → B
U + H → U
Next
ClrText
"(T, X, Y)"
{U, A, B} ◢
Menu "NEXT?", "YES", 0, "NO", 1
Lbl 1
"DONE"

Example

dx/dt = x sin y

dy/dt = y^2/500 + y/3 - y

Initial conditions: x(0) = 0.1, y(0) = 0.1, h = 0.1

T0 = 0
STEP = 0.1

ITERATIONS: 10
Results:  (T, X, Y): (1, 0.1075645239, 0.05134921355)

ITERATIONS: 10
Results:  (T, X, Y): (2, 0.1116714419, 0.02636554462)

Source:

John W. Harris and Horst Stocker.  Handbook of Mathematics and Computational Science.  Springer: New York.  2006.  ISBN 978-0-387-94746-4

Stack Exchange.  "Help with using the Runge-Kutta 4th order method on a system of 2 first order ODE'S"  Asked March 2014.   Last updated February 2019.  https://math.stackexchange.com/questions/721076/help-with-using-the-runge-kutta-4th-order-method-on-a-system-of-2-first-order-od  Retrieved June 17, 2020

Eddie

All original content copyright, © 2011-2020.  Edward Shore.   Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited.  This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author.

Sunday, July 26, 2020

HP 41C, HP 42S, TI-60: Arithmetic-Geometric Mean

HP 41C, HP 42S, TI-60:  Arithmetic-Geometric Mean



Arithmetic-Geometric Mean

The program AGM calculates the arithmetic-geometric mean of two positive integers x and y.   As the graphic above suggests, an iterative process is used to find the AGM, computing both the arithmetic mean and geometric mean until the two means converge.

a0 = x
g0 = y

Repeat:
Arithmetic Mean:  a1 = (a0 + g0)/2
Geometric Mean:  g1 = √(a0 * g0)
Transfer new to old:  a0 = a1, g0 = g1
Until |a1 - g1| < tolerance

You can set the tolerance as low as you want.  The programs presented on this blog set tolerance at 10^(-10)  (1E-10), to fit the calculator's display.

HP 41C Program: AGM

01 LBL^T AGM
02 STO 01
03 X<>Y
04 STO 02
05 X<>Y
06 LBL 00
07 RCL 02
08 RCL 01
09 ENTER
10 R↑
11 R↑
12 X<>Y
13 R↓
14 ENTER
15 R↑
16 +
17 2
18 /
19 STO 01
20 R↓
21 *
22 SQRT
23 STO 02
24 R↑
25 -
26 ABS
27 1E-10
28 X≤Y?
29 GTO 00
30 CLA
31 ^T AGM = 
32 ARCL 01
33 AVIEW
34 END

HP 42S/Swiss Micros DM42/Free42 Program AGM:

00 {53-Byte Prgm}
01 LBL "AGM"
02 STO 01
03 X<>Y
04 STO 02
05 X<>Y
06 LBL 00
07 RCL 02
08 RCL 01
09 ENTER
10 R↑
11 R↑
12 X<>Y
13 R↓
14 ENTER
15 R↑
16 +
17 2
18 /
19 STO 01
20 R↓
21 *
22 SQRT
23 STO 02
24 R↑
25 -
26 ABS
27 1E-10
28 X≤Y?
29 GTO 00
30 CLA
31 "AGM = "
32 ARCL 01
33 AVIEW
34 END

The instructions for both the HP 41C and 42S versions are same:  enter X and Y on the respective stacks and XEQ AGM.

Example (ALL/STD mode is applied):

AGM(37, 78): 
37, 78, XEQ AGM returns:
Alpha:  AGM = 55.5947005279

TI-60 Program: AGM

Instructions:

1.  Store X in memory register 1 and Y in memory register 2.
2.  Press [ RST ] [ R/S ], the value of |a1 - g1| is displayed.
3.  Keep on press [ R/S ] to repeat the calculation until |a1 - g1| falls under 10^(-10).
4.  Recall either memory register 1 or 2 to get the answer.

Registers needed: 1 - 4.

Step;  Key Code; Key
00;  71;  RCL
01;  01;  1
02;  85;  +
03;  71;  RCL
04;  02;  2
05;  95;  =
06;  55;  ÷
07;  02;  2
08;  95;  =
09;  61;  STO 
10;  03;  3
11;  71;  RCL
12;  01;  1
13;  86;  √
14;  65;  ×
15;  71;  RCL
16;  02;  2
17;  86;  √
18;  95;  =
19;  61;  STO
20;  04;  4
21;  71;  RCL
22;  03;  3
23;  61;  STO
24;  01;  1
25;  75;  -
26;  71;  RCL
27;  04;  4
28;  61;  STO
29;  02;  2
30;  95;  =
31;  87;  |X|
32;  13;  R/S

Example:

AGM(37, 78)
37 STO 1
78 STO 2
RST R/S

3.778495926, R/S
0.032100702, R/S
0.000002317, R/S
2 -11  (stop)

RCL 1 (or RCL 2):  55.59470053

Source:
"Arithmetic-geometric mean"  Wikipedia.  https://en.wikipedia.org/wiki/Arithmetic–geometric_mean  Last Edited June 12, 2020.  Accessed June 12, 2020.


Onward to August...

Eddie

All original content copyright, © 2011-2020.  Edward Shore.   Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited.  This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author.

Saturday, July 25, 2020

Fun with the 71B '20

Fun with the 71B '20

Differential Equations:  Runge Kutta Method 4th Order

Find a numerical solution to the differential equation:

dy/dx = f(x, y)

x is the independent variable, y is the dependent variable.  You define f(x,y) on line 10.   For example, dy/dx = sin(x * y) should have this as line 10:

10 DEF FNF(X,Y) = SIN(X*Y)

Line 5 is a remark line.  Remarks are followed by exclamation points on the HP 71B.

HP 71B Program: RK4 

Size: About 300 - 330 bytes

5 ! FNF(X,Y) = dY/dX
10 DEF FNF(X,Y) = [ insert f(x,y) here ]
15 DESTROY X,Y,H,K1,K2,K3,K4
20 INPUT "X0? "; X
25 INPUT "Y0? "; Y
30 INPUT "STEP? "; H
40 K1 = FNF(X,Y)
45 K2 = FNF(X+H/2,Y+H*K1/2)
50 K3 = FNF(X+H/2,Y+H*K2/2)
55 K4 = FNF(X+H,Y+H*K3)
60 X=X+H
65 Y=Y+H*(K1+2*K2+2*K3+K4)/6
80 DISP "(";X;",";Y;")" @ PAUSE
85 DISP "NEXT? Y/N"
90 A$=KEY$
95 IF A$="Y" THEN 40
99 IF A$="N" THEN DISP "DONE" @ END ELSE 85

Example:
dy/dx = sin(x * y) with inital condition y(0) = 0.5,  h = 0.2 * Ï€

First three results:
( .628318530718 , .607747199386 )   ( [ f ] [ +] (CONT), [ Y ] )
( 1.25663706144, 1.02432288082 )
( 1.88495559216, 1.51038862362 ) 

Hyperbolic Functions

[ S ] = sinh(x)
[ C ] = cosh(x)
[ A ] = asinh(x)
[ H ] = acosh(x)
[ X ]  to exit

HP 71B Program: HYP

Size:  401 bytes
acosh(x) requires that | x | ≥ 1

100 DESTROY A,X
115 DISP "sinh S/A, cosh C/H, X"
120 A$=KEY$
125 IF A$="S" THEN INPUT "X? ";X @ CALL SINH(X)
130 IF A$="C" THEN INPUT "X? ";X @ CALL COSH(X)
135 IF A$="A" THEN INPUT "X? ";X @ CALL ASINH(X)
140 IF A$="H" THEN INPUT "X? ";X @ CALL ACOSH(X)
145 IF A$="X" THEN 150 ELSE 115
150 DISP "DONE" @ END
200 SUB SINH(X)
205 DISP (EXP(X)-EXP(-X))/2 @ PAUSE
210 END SUB
300 SUB COSH(X)
305 DISP (EXP(X)+EXP(-X))/2 @ PAUSE
310 END SUB
400 SUB ASINH(X)
405 DISP LOG(X+SQR(X^2+1)) @ PAUSE
410 END SUB
500 SUB ACOSH(X)
510 DISP LOG(X+SQR(X^2-1)) @ PAUSE
515 END SUB

Example:
X = 2.86

sinh(2.86) returns 8.70212908815
cosh (2.86) returns 8.75939784845
asinh(2.86) returns 1.77321957441
acosh(2.86) returns  1.71190019325

Arithmetic-Geometric Mean

The arithmetic-geometric mean (AGM) is found by the iterative process:

a = 0.5 * (x + y)
g = √(x * y)

The values of a and g are stored into x and y, respectively.  The process repeats until the values of a and g converge.   A tolerance of 10^(-9) is used to display an 8-digit approximation.

HP 71B Program: AGM

Size:  148 Bytes

10 DESTROY X,Y,A,B
15 DISP "AGM(X,Y)" @ WAIT .5
20 INPUT "X? ";X
25 INPUT "Y? ";Y
30 A=.5*(X+Y)
35 G=SQR(X*Y)
40 X=A
45 Y=G
50 IF ABS(X-Y)>1E-9 THEN 30
55 DISP USING 60;X
60 IMAGE 10D.8D    // (10 digit integer parts with rounding to 8 decimal places)
65 END

Example:
AGM(178, 136)

Result:  156.29380544

Pythagorean Triple Generator

Given two positive integers m, n; where m > n, a Pythagorean triple is generated with the following calculations:

a = 2*m*n
b = m^2 - n^2
c = m^2 + n^2

Properties:

a^2 + b^2 = c^2
Perimeter: p = a + b + c
Area: r = a * b / 2

HP 71B Program: PYTHTRI

Size: 217 bytes

10 DESTROY M,N,A,B,C,R,P
20 DISP "M>N, INTEGERS" @ WAIT .5
25 INPUT "M? "; M
30 INPUT "N? "; N
35 A=2*M*N
40 B=M^2-N^2
45 C=M^2+N^2
50 P=A+B+C
55 R=A*B/2
60 DISP 'A = ';A @ PAUSE
65 DISP 'B = ';B @ PAUSE
70 DISP 'C = ';C @ PAUSE
75 DISP 'PERIM.=';P @ PAUSE
80 DISP 'AREA =';R
85 END

Example:
M = 16, N = 11

Results:
A = 352, B = 153, C = 377, P = 864, R = 23760

Impedance of An Alternating Current

The program ALTCURR calculates the impedance (magnitude and phase angle) of a sinusoidal alternating current consisting of one resistor, one capacitor, and one inductor in a series.

HP 71B Program: ALTCURR

Size: 210 bytes

10 DESTROY F,L,C, R,W,Z,T
15 DEGREES
20 INPUT "FREQUENCY? ";F
25 INPUT "INDUCTANCE? ";L
30 INPUT "CAPACITANCE? ";C
35 INPUT "RESISTANCE? ";R
40 W=2*PI*F
45 Z=SQR(R^2+(W*L-1/(W*C))^2)
50 T=ATAN((W*L-1/(W*C))/R)
55 DISP "MAGNITUDE= "; Z @ PAUSE
60 DISP "PHASE ANGLE = "; T

Example:
F = 152 Hz
L = 4.75E-3 H  (4.75 mH)
C = 8E-6 F  (8 μF)
R = 6400 Ω

Results: 
Magnitude:  6401.24704262
Phase Angle: -1.1309750812°

Source:
Rosenstein, Morton.  Computing With the Scientific Calculator Casio.  Japan. 1986.  ISBN 1124161430


Eddie

All original content copyright, © 2011-2020.  Edward Shore.   Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited.  This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author.

Sunday, July 19, 2020

Solving Systems of Quadratic Equations

Solving Systems of Quadratic Equations

This blog entry looks at four systems of quadratic equations.

In each system, A, B, C, D, E, and F are constants and we are trying to solve for x and y.  Note:  please verify each solution but substituting the computed x and y back in the equations. 

System 1

A* x^2 + B* y = C
D* x^2 + E *y = F

A * x^2 + B * y = C
-A/D * (D * x^2 + E * y) = F * -A/D

I'm going to solve for y first. 

A* x^2 + B * y = C
-A * x^2 - (A * E)/D * y = -(A * F)/D

(B - A * E / D) * y = C  - (A * F)/D

y = (C  - A * F / D) / (B - A * E / D)

Let y0 = y and now solve for x:

A * x^2 + B * y0 = C
A * x^2 = C - B * y0
x^2 = 1 / A * (C - B * y0)
x = ±√( 1 / A * (C - B * y0) )

Summary:

A* x^2 + B* y = C
D* x^2 + E *y = F

y = (C  - A * F / D) / (B - A * E / D)
x = ±√( 1 / A * (C - B * y) )

Example:

Pictures are generated by the HP Prime emulator.



x^2 + 5 * y = 10
4 * x^2 + 4 / 9 * y = 100

A = 1, B = 5, C = 10, D = 4, E = 4/9, F = 100

y = (C  - A * F / D) / (B - A * E / D)
y = (10 - 1 * 100 / 4) / (5 - 1 * 4/9 / 4)
y = -135 / 44 ≈ -3.068181818

x^2 = ±√( 1 / 1 * (10 - 5 * -135/44) )
x^2 = ±√( 1115 / 44 )
x ≈ ±5.033975476

The two points are:
( 5.033975476, -3.0681818)
( -5.033975476, -3.0681818)

System 2

A * x^2 + B * y^2 = C
D * x^2 + E * y^2 = F

We'll start with solving for x:

A * x^2 + B * y^2 = C
-B * D / E * x^2 - B * y^2 = -B * F / E

(A - B * D / E) * x^2 = C - B * F / E

x^2 = (C - B * F / E) / (A - B * D / E)

x = ±√( (C - B * F / E) / (A - B * D / E) )

Let x0 = x

A * x0^2 + B * y^2 = C
B * y^2 = C - A * x0^2
y^2 = 1 / B * (C - A * x0^2)

y = ±√( 1 / B * (C - A * x0^2) )

Summary:

A * x^2 + B * y^2 = C
D * x^2 + E * y^2 = F

x = ±√( (C - B * F / E) / (A - B * D / E) )
y = ±√( 1 / B * (C - A * x0^2) )

Example:



3 * x^2 + y^2 = 25/16
36 * x^2 + y^2 = 4

A = 3, B = 1, C = 25/16, D = 36, E = 1, F =4

x = ±√( (25/16 - 1 * 4 / 1) / (3 - 1 * 36 / 1) )
x = ±√( 13/176 ) ≈ ± 0.271778653

x^2 = 13/176

y = ±√( 1 / 1 * (25/16 - 3 * 13/176) )
y = ±√( 59/44 ) ≈ ± 1.157976291

Solutions:

(0.271778653, 1.157976291)
(0.271778653, -1.157976291)
(-0.271778653, 1.157976291)
(-0.271778653, -1.157976291)

System 3

A * x^2 + B * y^2 = C
D * x^2 + E * y = F

A * x^2 + B * y^2 = C
-A / D * (D * x^2 + E * y) = F * -A / D

A * x^2 + B * y^2 = C
-A * x^2 + -A * E / D * y = -A * F / D

B * y^2 + -A * E / D * y = C - A * F / D

B * y^2 + -A * E / D * y  + (A * F / D - C) = 0

Once y is solved for, then:

A * x^2 = C - B * y^2
x^2 = 1 / A * (C - B * y^2)
x = ±√ ( 1 / A * (C - B * y^2) )

Summary:

A * x^2 + B * y^2 = C
D * x^2 + E * y = F

B * y^2 + -A * E / D * y  + (A * F / D - C) = 0
x = ±√ ( 1 / A * (C - B * y^2) )

Example: 



49/16 * x^2 + 16 * y^2 = 64
x^2 + 7 * y = 5

A = 49/16, B = 16, C = 4, D = 1, E = 7, F = 5

B = 16
-A * E / D = -343/16
A * F / D - C = -779/16

y1 ≈ -1.19804377
x1 ≈ ±3.659362053

y2 ≈ 2.538548127
x2 ≈ ±3.573490854i

Solutions (real numbers):
( 3.659362053, -1.19804377)
( -3.659362053, -1.19804377)

System 4

A * x^2 + B * y^2 = C
D * x + E * y = F

E * y = F - D * x
y = F/E - D/E * x
y^2 = (D/E)^2 * x^2 - 2 * D * F / E^2 * x + (F/E)^2

A * x^2 + B * ((D/E)^2 * x^2 - 2 * D * F / E^2 * x + (F/E)^2) = C
A * x^2 + B * ((D/E)^2 * x^2 - 2 * D * F / E^2 * x + (F/E)^2) - C = 0
A * x^2 + B * (D/E)^2 * x^2 - (2 * B * D * F / E^2) * x + ( B * (F/E)^2 - C ) = 0

Once the solutions for x are found, solve for y by:
y = F/E - D/E * x

Summary:

A * x^2 + B * y^2 = C
D * x + E * y = F

(A + B * (D/E)^2) * x^2 + (-2 * B * D * F / E^2) * x + ( B * (F/E)^2 - C ) = 0
y = F/E - D/E * x

Example 4:




5 * x^2 + 6 * y^2 = 15
-10 * x + 13 * y = 2

A = 5, B = 6, C = 15, D = -10, E = 13, F = 2

A + B * (D/E)^2 = 1445/169
-2 * B * D * F / E^2 = 240/169
B * (F/E)^2 - C = -2511/169

x1 ≈ 1.2537792907
y1 ≈ 1.105994544

x2 ≈ -1.403882873
y2 ≈ -0.926063748

I hope you find this helpful.

Eddie

All original content copyright, © 2011-2020.  Edward Shore.   Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited.  This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author.

Python: Drawing a Ring of Circles

Python: Drawing a Ring of Circles The goal: Draw a number of circles that revolve around the origin (0,0). The center of each circle i...