Rationalizing a Quadratic Polynomial

Rationalizing a Quadratic Polynomial

Problem

Rewrite the Quadratic Polynomial 1 + A*x + B*x^2 as a rational function of polynomials.  A and B are real numbers (however, this should work if A and B were complex numbers).  

1 + A*x + B*x^2 → p(x) / q(x)

Attempt 1:  1 + A*x + B*x^2 → (1 + C*x) / (1 + D*x)

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

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

A*x + B*x^2 + D*x + A*D*x^2 + D*B*x^3 = C*x

Comparing the powers of x:

constant:  0 = 0

x:  A + D = C

x^2:  B + A*D = 0

x^3:  B*D = 0

This leads to either B=0 or D=0

Assume B=0. 

Then B + A*D = 0

A*D = 0

If A=0, then D = C, which leads to:

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

1 = (1 + C*x) / (1 + C*x)

1 = 1

If D=0, then A = C

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

1 + C*x = (1 + C*x) 

If D=0, then B=0, and we get the same results as above.   

Ultimately this transformation to (1 + C*x) / (1 + D*x) leads to nothing useful.

Attempt 2:  1 + A*x + B*x^2 → (1 + C*x^2) / (1 + D*x)

 

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

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

(A + D)*x + (B + A*D)*x^2 + B*D*x^3 = C*x^2

Comparing the powers of x:

constant:  0 = 0

x:  A + D = 0

x^2:  B + A*D = C

x^3:  B*D = 0

A + D = 0 implies that A = -D or D = -A

Also either B = 0 or D = 0.

Assume B = 0. Then with A = -D:

A*D = C

-D*D = C

C = -D^2

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

1 -  D*x = (1 - D*x^2) / (1 + D*x)

1 -  D*x = ((1 - D*x) * (1 + D*x))/ (1 + D*x)

1 - D*x = 1 - D*x

If we assume that D = -A, then:

C = -A^2 and

1 + A*x = (1 - A*x^2) / (1 - A*x)

1 + A*x = ((1 - A*x)  * (1 + A*x)) / (1 - A*x)

1 + A*x = 1 + A*x

Assume D = 0.

Then A = 0 and B = C:

1 + B*x^2 = 1 + B*x^2

1 + C*x^2 = 1 + C*x^2

Again, we have transformations that are trivial.

Attempt 3:  1 + A*x + B*x^2 → (1 + C*x^3) / (1 + D*x)

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

(1 + A*x + B*x^2) * (1 + D*x) = 1 + C*x^3

(A + D)*x + (B + A*D)*x^2 + B*D*x^3 = 1 + C*x^3

Comparing the powers of x:

constant:  0 = 0

x:  A + D = 0

x^2:  B + A*D = 0

x^3:  B*D = C

This implies that:

A + D = 0

D = -A

B + A*D = 0

B+ A*-A = 0

B = A^2   (this restricts A and B)

B * D = C

(A^2)*(-A) = C

C = -A^3

We can conclude that B = A^2, C = -A^3, D = -A

The relationship between A, B, C, and D are all connected in this case.

Examples:

A = 2 ⇒ B = 4, C = -8, D = -2 and

1 + 2*x + 4*x^2 = (1 - 8*x^3) / (1 - 2*x)

A = -2 ⇒ B = 4, C = 8, D = 2 and

1 - 2*x + 4*x^2 = (1 + 8*x^3) / (1 + 2*x)

Note:  Casio fx-991EX Week - September 5, 2022 to September 9, 2022 

Eddie

All original content copyright, © 2011-2022.  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. 

HP 17BII+ and TI-84 Plus CE Python: Zeta Approximation

HP 17BII+ and TI-84 Plus CE Python:  Zeta Approximation

Zeta Function

zeta(x) = Σ(1 ÷ (n^x), n =1 to ∞)

We are going to use the approximation:

zeta(x) = Σ(1 ÷ (n^x), n =1 to w) where w = intg(10^((a + 2) ÷ x)

where a is the number of decimal places desired.  The higher the accuracy, the longer the calculation takes.  Also the lower x is, the longer the calculation takes.

This is for all x > 0.

HP 17BII+ Formula ZETA

ZETA=0×L(W:10^IP((ACC+2)÷X))+Σ(N:1:G(W):1:INV(N^X))

Examples:

ACC =3, X = 2;  Result:  ZETA = 1.63

ACC =3, X = 3.5;  Result:  ZETA = 1.13

ACC =3, X = 8.7;  Result:  ZETA = 1.00

TI-84 Plus CE Python:  zeta.py

# 2021-12-07 ews

# zeta function approximation

from math import *

print("zeta function approximation")

x=eval(input("x? "))

a=eval(input("# places? "))

w=int(10**((a+2)/x)

z=0

n=1

while n<w:

  z+=(n**x)**-1

  n+=1

z=round(z,a)

print("zeta = "+str(z))

Eddie

All original content copyright, © 2011-2022.  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. 

Testing The Accuracy of LN and EXP Approximations

Testing The Accuracy of LN and EXP Approximations

Introduction

Today we are testing an approximation polynomial for two common scientific functions.  I used a TI-84 Plus CE in testing.   The values take a range from x = 0.5 to x = 11.  

Approximating e^x

Approximation used:

e^x ≈ 1 + x * (1 + x/2 * (1 + x/3 * (1 + x/4 * (1 + x/5 * (1 + x/6)))))

(fractions are simplified, see source)

TI-84 Plus CE Program:  EXPTEST

"EWS 2021-01-19"

Disp "e^(X) TEST","JON M. SMITH"

Disp "RANGE, MIN. 0.5"

Input "START? ",A

Input "INCREMENT? ",I

Input "NO STEPS? ",N

seq(X,X,A,A+I*N,I)→L₁

e^(L₁)→L₂

1+L₁*(1+L₁/2*(1+L₁/3*(1+L₁/4*(1+L₁/5*(1+L₁/6)))))→L₃

abs(L₂-L₃)→L₄

iPart(log(L₄))→L₅

Disp "L₁:X, L₂:EXP, L₃:APPROX","L₄:ABS ERROR","L₅:PLACES"

SetUpEditor L₁,L₂,L₃,L₄,L₅

Disp "PRESS STAT, 1. EDIT"

Test: x range:  0.5 to 11, increments of 0.5

Approximating ln x

Approximation used:

ln x ≈ y * (1 + y/2 * (1 + 2y/3 * (1 + 3y/4 * (1 + 4y/5))))

where x ≥ 0.5

and y = (x - 1)/x

(see source)

TI-84 Plus CE Program:  LNTEST

"EWS 2021-01-19"

Disp "ln(X) TEST","JON M. SMITH"

Disp "RANGE, MIN. 0.5"

Input "START? ",A

Input "INCREMENT? ",I

Input "NO STEPS? ",N

seq(X,X,A,A+I*N,I)→L₁

ln(L₁)→L₂

(L₁-1)/L₁→L₆

L₆*(1+L₆/2*(1+2*L₆/3*(1+3*L₆/4*(1+4*L₆/5))))→L₃

abs(L₂-L₃)→L₄

Disp "L₁:X, L₂:LN, L₃:APPROX","L₄:ABS ERROR"

SetUpEditor L₁,L₂,L₃,L₄

Disp "PRESS STAT, 1. EDIT"

Test: x range:  0.5 to 11, increments of 0.5

Final Thoughts

The approximation for ln x hold much better than e^x as x increases.   I would not recommend the following approximations for values above 5.  

Source:

Smith, Jon M.  Scientific Analysis on the Pocket Calculator John Wiley & Sons: New York. 1975

Eddie

All original content copyright, © 2011-2021.  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. 

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.

Casio fx-CG 50: Shortcuts Study

Casio fx-CG 50: Shortcuts Study

Introduction

When programming in calculators, sometimes it is useful to use alternate ways of accomplishing the task.  The alternate way could save keystrokes and memory, which can become important with older programming calculators and calculators where the number of steps available is about 150 or lower. 

Complex Number Functions

Some of the shortcut methods involves complex numbers and complex number operations, particular involving polar-rectangular conversions.

Below are the four common conversion calculations needed and alternate ways to accomplish them: 

√(x^2 + y^2) = abs(x + y ⅈ) = R>Pr(x,y) 

atan( y / x ) = arg(x + y ⅈ) = R>Pθ(x,y)**

r cos θ = real( r ∠ θ) = P>Rx(r,θ)

r sin θ = imag(r ∠ θ) = P>Ry(r,θ)

** The domain in the atan function is (-90°, 90°) while the arg (sometimes named atan2 and angle) and R>Pθ functions takes the (x,y) coordinates to account and the domain is (-180°, 180°).

Casio fx-CG50 Program: SHORTCUT

The program SHORTCUT illustrates six shortcut techniques regarding:

1.  Law of Cosines - Calculating the Third Side of a triangle

c = √(a^2 + b^2 - 2ab cos θ)
Complex Number Alternate: c = abs( real(a∠θ) - b + imag(a∠θ)*ⅈ )

2.  Rounding to the Nearest Integer

Alternate to the round function:  int( frac x + x 0

3.  √(1 - x^2),  |x| < 1

Trigonometric Alternate:  sin(acos x )

4.  x / √(1 + x^2),  any x

Trigonometric Alternate:  sin(atan x)

5.  √(a^2 + b^2 + c^2)

Complex Number Alternate:  abs(c + abs(a + b*ⅈ) * ⅈ)

6.  √(a^2 + 2*b^2)

Two Step Alternate (this is just one of many ways):
abs(a + b*ⅈ)
abs(b + ans*ⅈ)

Program:

Anything that follows // is a comment.  Don't type the comment in. 

a+bⅈ       // turn complex mode on
"BETTER PROG HP67/97"    // credit to source, slash character from CHAR menu
"EWS 2020-03-09"
Menu "SHORTCUTS","LAW OF COSINES",1,"ROUND TO INTEGER",2,
"√(1-X^2)",3," X÷√(1+X^2)",4,"√(A^2+B^2+C^2)",5,"√(A^2+2*B^2)",6

Lbl 1
"θ"?→θ
"A"?→A
"B"?→B
"√(A^2+B^2-2AB cos θ)"→Str 1
Exp(Str 1)→R            
"Abs(Rep(A∠θ)-B+Imp(A∠θ)ⅈ"→Str 2
Exp(Str 2)→S
Goto R

// Exp numerically evaluates a string 

Lbl 2
"NUMBER"?→X
"RndFix(X,0)"→Str 1
RndFix(X,0)→R      
// Exp(RndFix(a,b)) leads to an error message
"Int (Frac X+X)"→Str 2
Exp(Str 2)→S
Goto R

Lbl 3
"-1
"√(1-X^2) "→Str 1
Exp(Str 1)→R
"sin cos^-1 X"→Str 2
Exp(Str 2)→S
Goto R

Lbl 4
"X"?→X
"X÷√(1+X^2)"→Str 1
Exp(Str 1)→R
"sin tan^-1 X"→Str 2
Exp(Str 2)→S
Goto R

Lbl 5
"A"?→A
"B"?→B
"C"?→C
"√(A^2+B^2+C^2)"→Str 1
Exp(Str 1)→R
"Abs (C+Abs (B+A×ⅈ)×ⅈ)"→Str 2
Exp(Str 2)→S
Goto R

Lbl 6
"A"?→A
"B"?→B
"√(A^2+2B^2)"→Str 1
Exp(Str 1)→R
"Abs (A+Bⅈ):Abs (Ans+Bⅈ)"→Str 2
Exp(Str 2)→S
Goto R

Lbl R
ClrText
Black Locate(1,1,Str 1)
Black Locate(1,2,R)
Blue Locate(1,4,Str 2)
Blue Locate(1,5,S)

//  If you are using a monochrome screen such as the Casio fx-9750gII or fx-9860gII, leave out the Black and Blue commands.


You can download the program file here (fx-CG 50, 696 bytes):
https://drive.google.com/open?id=1b4PzPWEdQ0bV_l3l7vIcThVGxjOR94pb 

Source:
William Kob, John Kennedy, Richard Nelson.  Better Programming On The HP-67/97 PPC 1978

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.

HP 12C: Sums of x

HP 12C:  Sums of x

Introduction

The following are four programs to calculate:

*  ∑ x
*  ∑ x^2
*  ∑ x^3
*  ∑ ax + b

with limits x = 0 to n.

Each program uses a closed sum formula.  The second and fourth program uses Horner's Method. 

Note:  You can use the first three if your limits are x = 1 to n (because at x = 0, 0 is added to the sum).  However if you use the fourth program, to get the sum from 1 to n, subtract b. 

HP 12C Program:  ∑x (x=0 to n)

n
∑ x =  n * (n + 1) / 2
x=0

Enter n, execute program

STEP:  KEY CODE  KEY
01:  36  ENTER
02:  36  ENTER
03:  1  1
04:  40  +
05:  20  *
06:  2   2
07  10  ÷
08:  43,33,00  GTO 00

Keys Only:
ENTER
ENTER
1
+
2
*
÷
GTO 00

Example:  (n = 48)

48
∑ x  = 1176
x=0

HP 12C Program:  ∑x^2 (x=0 to n)

n
∑ x^2 =  n * ( ( n / 3 + 1 / 2 ) * n + 1 / 6 )
x=0

Enter n, execute program

STEP:  KEY CODE KEY
01:  36  ENTER
02:  36  ENTER
03:  36  ENTER
04:  3  3
05:  10  ÷
06:  2  2
07:  22  1/x
08:  40  +
09:  20  *
10:  6  6
11:  22  1/x
12:  40  +
13:  20  *
14:  43,33,00  GTO 00

Keys Only:
ENTER
ENTER
ENTER
3
÷
2
1/x
+
*
6
1/x
+
*
GTO 00

Example:  (n = 48)

48
∑ x^2  = 38,024
x=0

HP 12C Program:  ∑x^3 (x=0 to n)

n
∑ x^3 =  1 / 4 * (n^2 + n)^2
x=0

Enter n, execute the program

STEP:  KEY CODE  KEY
01:  36  ENTER
02:  36  ENTER
03:  2  2
04:  21  y^x
05:  40  +
06:  2:  2
07:  21  y^x
08:  4  4
09:  10  ÷
10:  43,33,00  GTO 00

Keys Only:
ENTER
ENTER
2
y^x
+
2
y^x
4
÷
GTO 00

Example:  (n = 32)

32
∑ x^3  = 278,784
x=0

HP 12C Program:  ∑(a*x + b)  (x=0 to n)

n
∑ (a*x + b) =  n * (a * n / 2 + (a / 2 + b) ) + b
x=0

Store a in R1, store b in R2, enter n on the x-stack, execute the program

STEP:  KEY CODE   KEY
01:  36   ENTER
02:  36   ENTER
03:  36   ENTER
04:  45, 1  RCL 1
05:  20   *
06:  2    2
07:  10   ÷
08:  45, 1  RCL 1
09:  2   2
10:  10   ÷
11:  45, 2   RCL 2
12:  40   +
13:  40   +
14:  20   *
15:  45, 2  RCL 2
16:  40  ÷
17:  43, 33, 00   GTO 00

Keys Only:
ENTER
ENTER
ENTER
RCL 1
*
2 
÷
RCL 1
2
÷
RCL 2
+
+
*
RCL 2
÷
GTO 00

Example:  (a = R1 = 3,  b = R2 = 4, n = 32)

3 STO 1,  4 STO 2,  32, R/S

32
∑ (3x + 4) = 1,716
x=0


I tried something new with my RPN keystroke programs, the second versions of each I just listed what keys to press without the step number or code.  I think this would be easier to read, but does leave out the key and key code.  Any comments on which style is preferred is appreciated. 


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.

HP Prime and HP 17BII: Trapezoid Rule Using Distinct Points

HP Prime and HP 17BII: Trapezoid Rule Using Distinct Points

Introduction

We can estimate the area of any surface by the use of sums and integral.  In calculus, we usually are given a function f(x), but here we are using measurements from one end to the other at various intervals. 

Technically, the intervals between each measurement do not have to be equal length.  However, having intervals of equal length makes things a lot easier, and in this blog entry, we assume they are. 

We have various methods to estimate the area.  One of the easiest ways is the Trapezoid Rule:

Area ≈ h/2 * ( y_1 + y_n + 2 * Σ( y_k , k, 2, n-1 ) )

Where:

h = interval length
y_k = length of each measurement, there are  n measurements
y_1 and y_n:  measurement of lengths at each end, respectively

Another rule to estimate area is the Simpson's Rule:

Area ≈ h/3 ( y_0 + y_n + 4 * Σ( y_k, k, 1, n-1, 2) + 2 * Σ( y_k, k, 2, n-2, 2) )

The program presented here uses the Trapezoid Rule. 

HP Prime Program AREAHGT

Two arguments:  h, a list of measurements

EXPORT AREAHGT(h, ms)
BEGIN
// h:  increment between measurements
// ms: list of measurements
// 2019-05-09 EWS
LOCAL k,n:=SIZE(ms);
RETURN h/2 * (ms(1) + ms(n) + 2 * Σ( ms(k), k, 2, n-1 ));
END;

HP 17BII+ (Silver)/HP 17BII Solver:  Trapezoid Rule

First:  define a SUM list named MS.  The solver uses that list to get the reference measurements.

Solver:

AREAHGT: AREA = 0 * L(N:SIZES(MS)) + H÷2 * (ITEM(MS:1) + ITEM(MS:G(N)) + 2 * Σ(K:2: G(N)-1: 1: ITEM(MS:K) )

Example

h = 0.5

MS:
y_1 = 1174
y_2 = 1078
y_3 = 979
y_4 = 984
y_5 = 810
y_6 = 779
y_7 = 800
y_8 = 852
y_9 = 966

Area:  3676

Eddie

All original content copyright, © 2011-2019.  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.

HP Prime and TI-84 Plus CE: Double Integral by Monte Carlo Method

HP Prime and TI-84 Plus CE:  Double Integral by Monte Carlo Method

Introduction

There are several ways to calculate integrals.  The most direct way is to find the anti derivative and evaluate it.  Another way is to use an approximation technique.   One technique to find a double integral is to evaluate random points on a given area (circle, rectangle, etc) and use the average in its calculations.

With the Monte Carlo method, we are at mercy of the random numbers selected.  Hence to improve our chances of getting a more accurate result, we have to use a large amount of random points for the calculation.  The flip side is that the Monte Carlo method, similar to the Romberg sum and Simpsons Method, finding the antiderivatve is not required. 

About MONTEINT

The program MONTEINT calculates the double integral

∫ ∫ f(x,y) dx dy

over a circular area with center (X, Y) and radius R.   Please be sure the function is continuous and exists over the entire circular region.   In order to help with the calculation time is that the program runs the calculations several times and averages the results.  You enter the number of batches (B, number of times the calculation is made) and the number of sample points of each batch (N). 

A first for me:  this is the first diagram I drew and posted using Microsoft Paint 3D! 

HP Prime Monte Carlo Double Integral Program: MONTEINT

EXPORT MONTEINT(f,C,D,R,B,N)
BEGIN
// 2018-12-11 EWS
// 'Z(X,Y)', center X,
// center Y, radius,
// number of batches,
// batch size

LOCAL V,S,J,K,I;

V:=0;
FOR K FROM 1 TO B DO
S:=0;
J:=0;
FOR I FROM 1 TO N DO
X:=RANDOM(C-R,C+R);
Y:=RANDOM(D-R,D+R);
IF (X-C)^2+(Y-D)^2≤R^2 THEN

// for EVAL() to work,
// X and Y must be global
S:=S+EVAL(f);
J:=J+1;
END;
END;
V:=V+(π*R^2/J)*S;
END;
V:=V/B;
RETURN V;
END;

Note:  Enter f(X,Y) in single quotes. 

TI 84 Plus Monte Carlo Double Integral Program: MONTEINT

Program MONTEINT

"2018-12-08 EWS"
Input "Z(X,Y) AS A STRING:",Y₀
FnOff 0
Disp "DISK"
Input "CENTER X: ",C
Input "CENTER Y: ",D
Input "RADIUS: ",R
"SETUP VOLUME"
0→V
Input "BATCHES: ",B
Input "SAMPLE SIZE: ",N
"LOOPS"
For(K,1,B)
"SUM, TRUE COUNT"
0→S
0→J
For(I,1,N)
2*R*rand+(C-R)→X
2*R*rand+(D-R)→Y
If (X-C)²+(Y-D)²≤R²
Then
S+Y₀→S
J+1→J
End
End
V+(π*R²/J)*S→V
Disp "BATCH",K,"/",B
End
V/B→V
Disp "APPROX DOUBLE","INTEGRAL:",V

Note:  I put code that shows progress since the TI-84 Plus CE has a slower processor than the HP Prime.

Examples

Full decimal answers will vary. 

Example 1:

f(X,Y) = √(LN(X + Y)) over a circle with center (5,5) of radius 1.   20 batches of 50 points each.

The value is approximately 4.76.

Example 2 - A Cautionary Tale:

f(X,Y) = X^2 + X*Y + Y^3 over a circle with center (0,0) of radius 2.  20 batches of 50 points each.

To illustrate how the Monte Carlo method is at the mercy of random points, the results vary between 11.3 to 13.1. 

I am going to increase the number of batches to 30 with 100 points each.  Thankfully with the HP Prime's processing speed, this won't increase the calculation by a large amount.  The TI-84 Plus CE will require a longer wait. 

After five trials, I get values between 12.07 and 13.15.

Increasing the number of points per batch to 200, I get results of 11.3 to 12.6.  The function f(X,Y) = X^2 + X*Y + Y^3  may not be a good function for the Monte Carlo method.

CAUTION:  The Monte Carlo is a method to calculate double integrals but it's only for approximate values only.  Several evaluations may be necessary to get a satisfactory approximation. 

Eddie

Source:

Ward Cheney and David Kincaid.  Numerical Mathematics and Computing Sixth Edition Thomson Brooks/Cole: Belmont, CA 2008  ISBN-13 978-0-495-11475-8

All original content copyright, © 2011-2018.  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.  Please contact the author if you have questions.

TI-84 Plus and HP Prime: Differential Equations and Half-Increment Solution, Numerical Methods

TI-84 Plus and HP Prime:  Differential Equations and Half-Increment Solution, Numerical Methods

Introduction

The program HALFSTEP solves the numerical differential equation

d^2y/dt^2 = f(dy/dt, y, t)  given the initial conditions y(t0) = y0 and dy/dt (t0) = dy0

In this notation, y is the independent variable and t is the dependent variable.

The Method

Let C = f(dy/dt, y, t).  Give the change of t as Δt.

First Step:

With t = t0:

h_1/2 = dy0 + C * Δt/2

y1 = y0 + dy0 * Δt

Loop:

t = t0 + Δt

h_I+1/2 = h_I-1/2 + C * Δt

y_I+1 = y_I +h_I+1/2 * Δt

Repeat as many steps as desired.

This method was presented by Robert M. Eisberg in his 1976 calculator programming book (see source below).

Variables

The program uses the following variables:

C:  d^2y/dt^2.   Represent dy/dt as the variable A, y as the variable Y, and t as the variable T.

The program will always designate Y as the independent variable and T as the dependent variable.

Examples:

Application	C	C for HALFSTEP
Free-Fall	d^2y/dt^2 = g	“9.80665” (SI) or “32.1740468” (US)
Free-Fall with Friction	d^2y/dt^2 = g - α (dy/dt)^2 (α = F/m)	“g - α * A^2” (sub numeric values for g, α)
Spring	d^2x/dt = -k/m * x	“-k/m * T” (sub numeric values for k, m)
Pendulum	d^2θ/dt = -α*sin(θ) (α = -g/l)	“-α * sin(Y)” (sub numeric values for α)
Damped, Driven Oscillations	d^2x/dt = -α*x – β*dx/dt + γ * sin(ω*t)	“-α*Y-β*A+γ*sin(ω*T)” (sub numeric values for α, β, γ)

HP Prime Program HALFSTEP

Input:  C.  Use single quotes to enclose d^2y/dt^2.  Represent dy/dt as A, y as Y, and t as T. 

Output:  A matrix of two columns, t and y.

EXPORT HALFSTEP(c,A,Y,D,tmax)

BEGIN

// d^2y/dt^2=C,dy0,y0,Δt,tmax

// EWS 2016-11-17

// C use single quotes

// 'dy=A, y=Y, t=T'

// Radian mode

HAngle:=0;

LOCAL mat:=[[0,Y]],T,H;

LOCAL K:=3,I;

T:=D;

H:=A+EVAL(c)*D/2;

Y:=Y+H*D;

mat:=ADDROW(mat,[D,Y],2);

FOR I FROM 2*D TO tmax STEP D DO

T:=I; A:=H;

H:=H+EVAL(c)*D;

Y:=Y+H*D;

mat:=ADDROW(mat,[I,Y],K);

K:=K+1;

END;

RETURN mat;

END;

TI-84 Plus Program HALFSTEP

Input:  For C, use enclose d^2y/dt^2 in quotes.  Represent dy/dt as A, y as Y, and t as T. 

Output:  A matrix of two columns, t and y.

"EWS 2016-11-27"

Func

Radian

Disp "D²Y/DT²=C"

Disp "USE A=DY/DT,Y,T"

Input "C, USE A STRING:",Y1

Input "DY0:",A

Input "Y0:",Y

Input "DELTA TIME:",D

Input "TIME MAX:",N

[[0][Y]]→[A]

D→T

A+Y1*D/2→H

Y+H*D→Y

augment([A],[[D][Y]])→[A]

For(I,2D,N,D)

I→T:H→A

H+Y1*D→H

Y+H*D→Y

augment([A],[[I][Y]])→[A]

End

[A]^T→[A]

Examples:

Please see the screen shots below.  Both are screen shots from the TI-84 Plus.

Source:  Eiseberg, Robert M.  Applied Mathematical Physics with Programmable Pocket Calculators  McGraw-Hill, Inc:  New York.  1976.  ISBN 0-07-019109-3

This blog is property of Edward Shore, 2016.

HP Prime and HP 49G/50G: Tridiagonal Matrices

HP Prime and HP 49G/50G:  Tridiagonal Matrices

The program TRIDIAG will build a tridiagonal matrix of the form:

D1	R1	0	0	…	0	0	0
L1	D2	R2	0	…	0	0	0
0	L2	D3	R3	…	0	0	0
…	…	…	…	…	…	...	…
0	0	0	0	…	L(n-2)	D(n-1)	R(n-1)
0	0	0	0	…	0	L(n-1)	D(n)

Where

D represents the list/vector of diagonal elements, and D has the length n,

L represents the list/vector of elements to the left of the diagonal, and L has the length n-1, and

R represents the list/vector of elements to the right of the diagonal, and R has the length of n-1.

HP Prime Program tridiag

Input:   tridiag(ll, ld, lr) where the arguments are lists.

Program:

EXPORT tridiag(ll, ld, lr)

BEGIN

LOCAL m,s,j,t;

s≔size(ld);

m≔diag(ld);

FOR j FROM 1 to s-1 DO

t≔j+1;

m(t,j)≔ll(j);

m(j,t)≔lr(j);

END;

RETURN m;

END;

HP 49G/50g Program TRIDIAG

This is the first program I made on a HP 49G, but it should work on the 50g just fine.

Input:

3:  vector of left elements

2:  vector of diagonals

1:  vector of right elements

Program:

≪ → L D R

≪ D SIZE OBJ→

DROP DUP D

SWAP DIAG→ →

S M

≪ M 1 S 1 –

FOR J

J 1 + J 2 →LIST

L J GET PUT

J J 1 + 2 →LIST

R J GET PUT

NEXT ‘J’ PURGE ≫ ≫ ≫

Example:

Left elements:  4, -3, 8

Diagonals: 1, 2, 3, 4

Right elements: -2, 5, 1

There are 4 diagonal elements, therefore the matrix will be 4 x 4.

Inputs:

HP Prime:  tridiag({4,-3,8},{1,2,3,4},{-2,5,1})

HP 49G/50g (RPN Mode):

3:  [4, -3, 8]

2:  [1, 2, 3, 4]

1:  [-2, 5, 1]

Result:

1	-2	0	0
4	2	5	0
0	-3	3	1
0	0	8	4

I hope you like my new style for matrices. 

Happy Halloween!

This blog is property of Edward Shore, 2016.