Numworks: Allowing Repeated Calculations in Python

Numworks: Allowing Repeated Calculations in Python

Introduction

Say we want the user to repeat a calculation or a routine for as long as the user wants. At the end of each calculation, the user is asked if they want another calculation. The user enters one of two possibilities:

1 for an additional calculation

0 for no more calculation

A flag variable is set up with a default value of 1. As long as this value remains at 1, the calculation repeats.

Scrip structure (I named this variable rptflag, but it can be any name):

rptflag=1

while rptflag!=0:

  …

  < insert inputs and show results >

  …

  print(“Again? No = 0, Yes = 1”)

  rptflag=eval(input())

Entering 0 for the Again? Question exits the loop. We need the eval (evaluate) or int (input) because leaving the input command alone defaults the input to a string.

For the sample scripts below, the escape velocity off the surface of planets and other celestial objects are calculated.

Repeated Loop: General

This script does not require any special modules.

Script: rptdemo1 (We can use this on any calculator or platform)

from math import *

# repeat demonstration 

# set up repeat flag

rptflag=1

while rptflag!=0:

  print("ESCAPE VELOCITY")

  m=eval(input("Mass (kg)? "))

  r=eval(input("Radius (m)? "))

  v=sqrt(1.33486e-10*m/r)

  print("Velocity: ",v,"\nm/s")

  print("Again? no=0, 1=yes")

  rptflag=eval(input())

Repeated Loop: Numworks – using the keyboard

Instead of having the user enter 0 or 1, the user presses the keys 0 or 1. To allow the user key presses, use the Ion module. KEY_ZERO is for the 0 key and KEY_ONE is for the 1 key.

To make the screen look presentable, I’m clearing the screen. This is accomplished by the use of the Kandinsky module. The use of the Kandinsky module was suggest to me on Reddit by Feeling_Walrus3033, and I give credit and gratitude for the suggestion.

The line:

fill_rect(0,0,320,240,(255,255,255)) creates a blank white screen.

Clearing the screen this way will call for a way to display text on the screen. Unfortunately the print() command won’t do it. Kandinsky comes to the rescue with the draw_string command. The syntax for the draw_string command is:

draw_string(“text” or variable containing the string, x pixel, y pixel, [text color], [background color])

Colors are optional. The default is black text on white background.

Keep in mind that using this method will require more memory.

Script: rptdemo2 (Specific to Numworks, for other platforms and calculator, check your specific manual)

from math import *

from ion import *

from kandinsky import *

# repeat demonstration 

# set up keys

def key():

  while True:

    if keydown(KEY_ZERO):

      return 0

    if keydown(KEY_ONE):

      return 1

# set up repeat flag

rptflag=1

while rptflag==1:

  print("ESCAPE VELOCITY")

  m=eval(input("Mass (kg)? "))

  r=eval(input("Radius (m)? "))

  v=sqrt(1.33486e-10*m/r)

  

  # build strings

  s1="Velocity: "+str(v)+" m/s"

  s2="Again? no=0, 1=yes"

  fill_rect(0,0,320,240,(255,255,255))

  draw_string(s1,0,0)

  draw_string(s2,0,20)

  rptflag=key()

draw_string("DONE",0,40,(255,0,0))

Sample Data to Try

	Mass (kg)	Radius (m)	Escape Velocity (m/s)
Earth	5.972168E24	6378.137E3	11179.88618486758
Jupiter	1.8982E27	71492E3	59533.32250562
Sun (Solar System)	1.9885E30	6.957E8	617688.6988877566

E: [ ×10^x ] button, shown as “e”

Note: I will be in Nashville for the 2024 HP Handheld Conference on September 21-22, 2024. The next scheduled post will be on September 28, 2024.

Take care,

Eddie

All original content copyright, © 2011-2024. 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 CAS: Riemann-Louiville Integral vs Taking the Indefinite Integral Twice

HP Prime CAS: Riemann-Louiville Integral vs Taking the Indefinite Integral Twice

Introduction

The Riemann-Louiville Integral takes the integral of the function f(x) of any positive order v. The integral is defined as:

c_D_x^(-v) = 1 / Γ(v) * ∫( (x – t) * f(t) dt, t = c, t = x)

where:

c = a real constant, which can be zero

f(x) = function of x

t = dummy variable of integration

‘v = order where v >0

If v=1, this is the regular integral. However, the value of v can be a positive non-integer. If v=2, then the Riemann-Louiville integral is a result if you integrate the function twice.

This is one of the formulas on determining indefinite integrals of various orders.

HP Prime CAS Function: dblint

Double Integral of f(x) which takes the integral of f(x) twice. The variable x is used in the function.

dblint(f):= ∫∫ f dx dx

HP Prime CAS Function: rlint

The Riemann-Liouville Integral of f(x). The input has the variable x as the independent variable. The result of the function returns t as the independent variable.

rlint(f,c,v):=(∫(t–x)^(v–1)*f,x,c,t)) / Gamma(v)

Note that the variables x and t are switched to allow the input to be a function of x.

Notes

This was programmed on the CAS page in the format:

func(var) := function

I was not able to use the Program Editor mode at time of programming (July 30, 2024). (Beta Firmware 15048)

Examples

	Double Integral: dblint	RLI, v = 2: rlint with c = 0
f(x) = x^m, m>0	x^(m+2) / (m^2 + 3*m + 2)	t^(m+2) / (t^2 + 3*t + 2)
f(x) = a*x + b	(a*x^3 + 3*b*x^2) / 6	(a*t^3 + 3*b*t^2) / 6
f(x) = e^x	e^x	-t + e^t – 1
f(x) = cos x	-cos x	-cos t + 1

Taking the derivative twice will return us back to the original function. Note that the indefinite integral function assumes that the added constant is zero ( ∫ f(x) dx = F(x) + C ).

Source

Kimeu, Joseph M., "Fractional Calculus: Definitions and Applications" (2009).Masters Theses & Specialist Projects. Paper 115. http://digitalcommons.wku.edu/theses/115

Eddie

All original content copyright, © 2011-2024. 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 Numworks: Plotting a Parametric Line of Motion

 HP Prime and Numworks: Plotting a Parametric Line of Motion

Plotting the Position of Motion

This program draws a 2D motion plot from an initial starting point (x, y) given initial velocity and acceleration. The rate and direction of both velocity and acceleration are assumed to be constant.

x(t) = ax * t^2 / 2 + vx * t + x0

y(t) = ay * t^2 / 2 + vy * t + y0

where:

t = number of seconds.

ax = acceleration in the x direction

ay = acceleration in the y direction

vx = initial velocity in the x direction

vy = initial velocity in the y direction

The HP Prime PPL program PLOTMOTION displays a traceable curve with a table of values.

The Numworks python script plotmtn.py uses math and matplot.pyplot modules. A scatter plot is laid on top of the path plot. The plot begins at the initial point where it is marked green.

HP Prime Program: PLOTMOTION

EXPORT PLOTMOTION()

BEGIN

// EWS 2024-07-22

LOCAL ch;

ch:=INPUT({C,D,A,B,V,U,N},

"Motion Plot Per Second",

{"x0:","y0:","vx:","vy:","ax:","ay:","n:"},

{"initial x position",

"initial y position",

"initial velocity x direction",

"initial velocity y direction",

"acceleration x direction",

"acceleration y direction",

"number of seconds"});

// user presses cancel

IF ch==0 THEN

KILL;

END;

// user presses OK

STARTAPP("Parametric");

'V*T^2/2+A*T+C'▶X1;

'U*T^2/2+B*T+D'▶Y1;

CHECK(1);

Parametric.Tmin:=0;

Parametric.Tmax:=N;

Parametric.Tstep:=1;

// table and plot

STARTVIEW(10);

STARTVIEW(9);

END;

Numworks Python Code: plotmtn.py

from math import *
from matplotlib.pyplot import *

# 2024-07-23 EWS

print("Motion Plot per second from (0,0)")
c=eval(input("init. x position? "))
d=eval(input("init. y position? "))
a=eval(input("init. x velocity? "))
b=eval(input("init. y velocity? "))
v=eval(input("acceleration x? "))
u=eval(input("acceleration y? "))
n=int(input("number of seconds? "))

x=[v*t**2/2+a*t+c for t in range(n)]
y=[u*t**2/2+b*t+d for t in range(n)]

x0=min(x)-1
x1=max(x)+1
y0=min(y)-1
y1=max(y)+1
axis((x0,x1,y0,y1))
text(x0,y1-1,"Motion Plot")
plot(x,y,'gray')
scatter(x,y,color='blue',marker="h")
scatter(x[0],y[0],color='green',marker="h")
show()

Source

Tremblay, Christopher. Mathematics for Game Developers. Thomson Course Technology. Boston, MA. 2004. ISBN 1-59200-038-X.

Eddie

All original content copyright, © 2011-2024. 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: Minimum Distance Between a Point and a Line

HP Prime: Minimum Distance Between a Point and a Line

Introduction

We have a line in the form of y = m * x + b, where m is the slope of the line and b is the y-intercept of the line, and a separate point (px, py). The task is to find the minimum distance, or the shortest distance, between the point and the line. The separate point is not required to be on the line. The line and point are in two-dimensional space.

If the point (px, py) is not on the line, then theoretically, there are an infinite amount of distances between the point and the line. However, to get the shortest distance, draw a path that is “directly straight” to the line. This is achieved by choosing a line that connects the (px, py) that is a line that is orthogonal (perpendicular) to the line y = m * x + b.

The line drawn is of the form y = -1/m * x + b1. The slope of the orthogonal line is -1/m. Assuming that m ≠ 0, the y-intercept of the orthogonal line is b1 = y1 + x1 / m.

The next step is to find where the two lines intersect, which is done by solving the following system for x and y:

y = m * x + b

y = -m / n + b1

Label the intersection point (x1, y1). The minimum distance will be calculated as follows:

dist = √( (x1 – px)^2 + (y1 – py)^2 ) = abs( (x1 – px) + (y1 – py)*i)

If m = 0, the line is in the form of y = b. The orthogonal line is x = px, and the distance is simply abs( (y1 – py)*i ).

HP Prime Code: PTLNDIST

EXPORT PTLNDIST()

BEGIN

// 2024-07-21 EWS

// radian

HAngle:=0;

LOCAL px,py,m,b;

INPUT({m,b,px,py},

"Point-Line Distance (px, py), y=mx+b",

{"m:","b:","px:","py:"},

{"m: slope"," b: y-intercept",

"point x","point y"});

LOCAL y0;

y0:=m*px+b;

LOCAL b1,mt,x1,y1,dist,str;

IF m≠0 THEN

b1:=py+px/m;

mt:=[[−m,1],[1/m,1]]^-1*[[b],[b1]];

x1:=mt[1,1];

y1:=mt[2,1];

dist:=ABS((x1-px)+(y1-py)*√(-1));

ELSE

x1:=px;

y1:=b;

dist:=ABS((y1-py)*√(-1));

END;

// print results

PRINT();

PRINT("Results:");

PRINT("Intersect point:");

PRINT("x: "+STRING(x1));

PRINT("y: "+STRING(y1));

PRINT("");

IF m≠0 THEN

str:="Y="+STRING(-1/m)+"*X+"+STRING(b1);

ELSE

str:="X="+STRING(x1);

END;

PRINT("Orthogonal Line:");

PRINT(str);

PRINT("");

PRINT("Minimum Distance:");

PRINT(dist);

RETURN {x1,y1,str,dist};

END;

Note:

√(-1) represents the imaginary number ⅈ ( [ Shift ], [ 2 ] ).

Inputs:

* The slope of the y-intercept of the line y = m * x + b (no vertical lines, but m can be zero)

* The point (px, py)

Outputs:

* The line that runs through point (px, yx) that is orthogonal to y = m * x + b. The slope and y-intercept of the orthogonal line, which the line will be stated in a string

* The intersection point of the two lines.

* The distance between (px, yx) and the intersection point. (dist)

Examples

Example 1:

Inputs: Line: y = 5 x – 2, Point: (-1, -5)

m = 5

b = -2

px = -1

py = -5

Results:

Intersect point:

x = -0.615384615386

y = -5.07692307692

Orthogonal Line:

Y = -0.2 * X – 5.2

Minimum distance:

0.392232270274

Example 2:

Inputs: Line: y = 6, Point: (3, -9)

m = 0

b = 6

px = 3

py = -9

Results:

Intersect point:

x = 0.764705882353

y = 4.05882352941

Orthogonal Line:

X = 3

Minimum distance:

15

Example 3:

Inputs: Line: y = 4 x + 1, Point: (5, 3)

m = 4

b = 1

px = 5

py = 3

Results:

Intersect point:

x = 0.764705882353

y = 4.05882352941

Orthogonal Line:

Y = -0.25 * X + 4.25

Minimum distance:

4.36564125066

Source

Tremblay, Christopher. Mathematics for Game Developers. Thomson Course Technology. Boston, MA. 2004. ISBN 1-59200-038-X.

Eddie

All original content copyright, © 2011-2024. 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 Update: Official Firmware 15157

 HP has released an official firmware for the HP Prime:  Firmware 15157.  When you start the HP Connectivity Kit, the program should prompt you to download the firmware.

More information and downloads from TI-Planet (the page is in French):

https://tiplanet.org/forum/viewtopic.php?p=276332#p276332

Eddie

All original content copyright, © 2011-2024. 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 and Python: Gaussian Quadrature

 Casio fx-CG 50 and Python:  Gaussian Quadrature

The Method

Gaussian quadrature is a method of numerical integration of a function f(x) over an interval [a, b] where a < b. Quadrature works similar to the Trapezoid Rule or Simpson's Rule, but instead of equally spaced points of measure, points are determined by the roots of a certain polynomial. Polynomials used include the Legendre, Jacobi, Chebyshev, Hermite, and Laguerre. Each point, xi, will be assigned a corresponding certain weight, wi, designed to make the approximation the most accurate as possible.

In general, ∫( f(x) dx, x = a to x = b) ≈ Σ( wi * f(xi) from i = 1 to n)

For today’s blog, the Legendre polynomial is used and the domain of integration is restricted to the interval [ -1, 1 ]. Hence:

∫( f(x) dx, x = -1 to x = 1) ≈ Σ( wi * f(xi) from i = 1 to n)

Steps:

1. Determine the number of points (order) is needed.

2. Find the roots of the Legendre polynomial, P_n(x) = 0. This determines the points xi.

3. Calculate the weights as follows:

wi = 2 / ((1 – xi^2) * (P’_n(xi))^2)

where P’_n(xi) is the numerical derivative of the Legendre polynomial at point x = xi.

4. Approximate ∫( f(x) dx, x = -1 to x = 1) ≈ Σ( wi * f(xi) from i = 1 to n).

Legendre Polynomials

There are several ways to determine the Legendre polynomials of order n, these are just a few:

P_n(x) = 1 / (2^n * n!) d^n/dx^n (x^2 – 1)^n

P_n(x) = Σ( comb(n, k), * comb(n + k, k) * comb(x -1, 2)^k for k = 0 to n)

where comb is the combination method

Recursive method:

(n + 1) * P_n(x) = (2*n + 1) * x * P_n-1(x) – n * P_n-2(x)

where P_0(x) = 1 and P_1(x) = x.

A table of Legendre polynomials of orders 0 to 10 is found here (scroll down a bit):

https://en.wikipedia.org/wiki/Legendre_polynomials

An Example: Order 3

Let’s determine the points xi and the weights wi for order 3 (n = 3).

The Legendre polynomial of order 3 is:

P_3(x) = 1 / 2 * (5 * x^3 – 3 * x) = 5 / 2 * x^3 – 3 / 2 * x

And it’s derivative:

P’_3(x) = 15 / 2 * x^2 – 3 / 2

Solve P_3(x) = 0 yields:

5 / 2 * x^3 – 3 / 2 * x = 0

x * (5 / 2 * x^2 – 3 / 2) = 0

x= 0 or 5 / 2 * x^2 – 3 / 2 = 0:

5 / 2 * x^2 – 3 / 2 = 0

5 / 2 * x^2 = 3 / 2

x^2 = 3 / 5

x = ±√(3 / 5)

The roots are, from lowest value to highest value are: x1 = - √(3 / 5), x2 = 0, x3 = √(3 / 5)

Determining the associated weights yields: w1 = 5 / 9 , w2 = 8 / 9, w3 = 5 / 9

Our points and weights are:

x1 = -√(3 / 5), w1 = 5 / 9

x2 = 0, w2 = 8 / 9

x3 = √(3 / 5), w3 = 5 / 9

Casio fx-CG 50 Program: QUADRAT0

(The last character in the file name is a zero)

“3PT QUADRATURE”

“-1 TO 1”

“NO QUOTES”

“F(X)”? → fn1

∫(fn1, A, B) → D

{-√(3 _| 5), 0, √(3 _| 5) } → List 1

{ 5 _| 9, 8 _| 9, 5 _| 9 } → List 2

0 → S

For 1 → I To 3

List 1[ I ] → X

List 2[ I ] → W

S + W × fn1 → S

Next

ClrText

Black Locate 1, 3, “∫(fn1 DX)”

Blue Locate 1, 4, “D=”

Red Locate 4, 4, D

Blue Locate 1, 5, “S=”

Red Locate 4, 5, S

_| is the fraction symbol, pressed by the [ [] / [] ] key.

If you have a monochrome calculator, such as the fx-9750G series or the fx-9860G series, leave out the color commands (Black, Blue, Red) in front of the Locate commands.

The function is stored in function memory, slot 1, which is found in the OPTN menu.

Variables:

D = actual integral, using the integral command of the calculators

S = approximation determined the quadrature

Examples

Radians angle assumed	D (actual)	S (approximate)
∫( sin x dx, x = -1 to x = 1)	0	0
∫( e^(.5 * x) dx, x = -1 to x = 1)	2.084381222	2.084380222
∫( ln(x + 4) / 2 dx, x = -1 to x – 1)	1.365058335	1.365060354

7^th Order and Expanding the Range

The 7^th order Legendre polynomial is:

P_7 = 1 / 16 * (429 * x^7 – 693 * x^5 + 315 * x^3 – 35 * x)

The associated roots (xi) and weights (wi) are:

Point, xi =	Weight, wi =
-0.949107912342759	0.129484966168870
-0.741531185599394	0.279705391489277
-0.405845151377397	0.381830050505119
0	0.417959183673469
0.405845151377397	0.381830050505119
0.741531185599394	0.279705391489277
0.949107912342759	0.129484966168870

These values are from the table presented in “Gauss-Kronrod quadrature formula” (see Source). I verified the roots with Wolfram Alpha, which you can see here:

https://www.wolframalpha.com/input?i2d=true&i=Divide%5B1%2C16%5D+*++%5C%2840%29429+*+Power%5Bx%2C7%5D+%E2%80%93+693+*+Power%5Bx%2C5%5D+%2B+315+*+Power%5Bx%2C3%5D+%E2%80%93+35+*+x%5C%2841%29%3D0

We are not limited to the interval [-1, 1]. To use any interval [a, b], we can scale the points as so:

xi ‘ = (b – a) / 2 * xi + (b + a) / 2

wi’ = (b – a) / 2 * wi

Casio fx-CG50 Program: GAUSS7

This the 7-point Gaussian quadrature, the calculator basic version.

“7PT GAUSS RULE”

“NO QUOTES”

“F(X)”? → fn1

“A”? → A

“B”? → B

∫(fn1, A, B) → D

{ -0.949107912342759,-0.741531185599394,-0.405845151377397,0,0.405845151377397,0.741531185599394,0.949107912342759 } → List 1

{ 0.12948496618870,0.279705391489277,0.381830050505119,0.417959183673469,0.381830050505119,0.279705391489277,0.129484966168870 } → List 2

0 → S

For 1 → I To 7

(B – A) ÷ 2 × List 1[ I ] + (B + A) ÷ 2 → X

List 2[ I ] × (B – A) ÷ 2 × fn1 + S → S

Next

ClrText

Black Locate 1, 3, “∫(fn1 DX)”

Blue Locate 1, 4, “D=”

Red Locate 4, 4, D

Blue Locate 1, 5, “S=”

Red Locate 4, 5, S

Python Code: gauss7.py

The only module needed is the math module, so it should work with all calculators. This script was created with the fx-CG 50.

from math import *

# 7 point Gauss Integral Approximation

def f(x):

return (x-1/2)*(x+1/4)

# xi and wi

xi=[-0.949107912342759,-0.741531185599394,-0.405845151377397,0,0.405845151377397,0.741531185599394,0.949107912342759]

wi=[0.12948496618870,0.279705391489277,0.381830050505119,0.417959183673469,0.381830050505119,0.279705391489277,0.129484966168870]

# sum

s=0

# limits

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

b=eval(input("b? "))

# integral approximation

for i in range(7):

  x=(b-a)/2*xi[i]+(b+a)/2

  w=(b-a)/2*wi[i]

  s+=w*f(x)

print("Approx. Integral:\n",s)

Examples

Radians mode	D (actual)	S (Casio Basic)	S (Python)
∫( 1 / √(x + 3) dx, 1, 5)	1.656854249	1.658042205	1.656854249504868
∫( .03 * e^(-x^2) dx, -1, 1)	0.04480944797	0.04484387349	0.04480944866632411
∫( e^(-x) dx, -5, 0)	147.4131591	147.5268171	147.4131590824608
∫( .3 * ( cos x – 0.05) dx, 0, π / 2)	0.2764380551	0.2767068146	0.2764380551025117
∫( (x – 1 / 2) * (x + 1 / 4) dx, -1, 1)	0.4166666667	0.4168576653	0.41666666666867568

Notice a difference between the Basic and Python?

Observations

* The program calls for approximations of xi and wi that are carried out 15 digits. Since the roots of the Legendre polynomials are most likely to be irrational, the more digits used, the better (generally).

* Most calculators will store up to 10 to 13 digits internally.

* By contrast, Python will store many more decimal places, at least 15 to 17.

* The approximation uses a lot of calculations. Round off errors can affect the final results, even when the round off is small.

* I recommend the Quadrature method if you are using Python or any other platform that allows for a lot of decimal places.

* Thankfully, most calculators have the numeric integration function.

Sources

“Guass-Kronrod quadrature formula” Wikipedia. Last Edited December 27, 2023. Accessed July 3, 2024. https://en.wikipedia.org/wiki/Gauss%E2%80%93Kronrod_quadrature_formula

“Gaussian quadrature” Wikipedia. Last Edited July 7, 2024. Accessed July 7, 2024. https://en.wikipedia.org/wiki/Gaussian_quadrature

This is a long one, folks. Thank you and have a great day. To all the students, may your year of studies be full of success, discovery, and joy.

Until next time,

Eddie

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