RPN: HP 32S/DM32: Simple Harmonic Motion

RPN: HP 32S/DM32: Simple Harmonic Motion

Introduction

The simple harmonic motion (SHM) of a spring can be described by Hooke’s Law:

F = -K * x(t)

x(t) = position of the spring

F = force of the spring

Note that: F = mass * acceleration = M * x’’(t)

Then:

M * x’’(t) = - K * x(t)

M * x’’(t) + K * x(t) = 0

x’’(t) + (K / M) * x(t) = 0

Without going into much details, two general solutions can have the form:

Case I:

x(t) = A * sin(W * t + B) with W = √(K / M) (A, B are constants)

and

Case II:

x(t) = A * cos(W*T + D) with W = √(K / M) (A, D are constants)

Believe it or not, both solutions work.

Case I:

x(t) = A * sin(W * t + B)

x’(t) = A * W * cos(W * t + B)

x’’(t) = -A * W^2 * sin(W * t + B)

Observe:

x’’(t) = -A * W^2 * sin(W * t + B)

x’’(t) = -W^2 * (A * sin(W * t + B))

x’’(t) = -W^2 * x(t)

W = √(K / M), W^2 = K / M

x’’(t) = -(K / M) * x(t)

M * x’’(t) = -K * x(t)

We are back to Hooke’s Law.

Case II:

x(t) = A * cos(W * t + D)

x’(t) = -A * W * sin(W * t + D)

x’’(t) = -A * W^2 * cos(W * t + D)

Then:

x’’(t) = -W^2 * (A * cos(W * t + D))

x’’(t) = -W^2 * x(t)

x’’(t) = -(K / M) * x(t)

M * x’’(t) = -K * x(t)

We once again arrive at Hooke’s Law. Another valid solution.

These two solutions are chosen for simplicity. A general solution would take the form of:

x(t) = A2 * cos(W * t + D) + A2 * sin(W * t + B)

Determining the Constants

Let:

x0 = initial position, at t = 0

v0 = initial velocity, at t = 0

Let’s use the simplified solutions.

Case I: Sine Approach

x(t) = A * sin(W * t + B)

x’(t) = A * W * cos(W * t + B)

x0 = A * sin(B)

v0 = A * W * cos(B)

x0 = A * sin(B)

v0 / W = A * cos(B)

Square each equation and add:

x0^2 + (v0 / W)^2 = A^2 * ((sin B)^2 + (cos B)^2)

x0^2 + (v0 / W)^2 = A^2 * 1

√( x0^2 + (v0 / W)^2 ) = A

A = √( x0^2 + (v0 / W)^2 )

Divide the top equation by the bottom equation:

x0 / (v0 / W) = tan(B) (assume A ≠ 0)

B = arctan( x0 / (v0 / W) )

⇒ atan2( x0, (v0 / W)) = arg( (v0 / W) + x0*i ) where i = √-1

Case II: Cosine Approach

x(t) = A * cos(W * t + D)

x’(t) = -A * W * sin(W * t + D)

x0 = A * cos(D)

v0 = -A * W * sin(D)

x0 = A * cos(D)

v0 / W = -A * sin(D)

Square each equation and add:

x0^2 + (v0 / W)^2 = (A * cos(D))^2 + (-A * sin(D))^2

x0^2 + (v0 / W)^2 = A^2 * cos(D)^2 + (-A) *(-A) * sin(D)^2

x0^2 + (v0 / W)^2 = A^2 * cos(D)^2 + A^2 * sin(D)^2

x0^2 + (v0 / W)^2 = A^2 * 1

√( x0^2 + (v0 / W)^2 ) = A

A = √( x0^2 + (v0 / W)^2 )

Divide the top equation by the bottom equation:

x0 / (v0 / W) = -cot(D) (assume A ≠ 0)

(v0 / W) / x0 = -tan(D)

-(v0 / W) / x0 = tan(D)

D = arctan(-(v0 / W) / x0)

⇒ atan2( -(v0 / W), x0) = arg( x0 - (v0 / W)*i)

The atan2, arg, and similarly, the rectangular to polar conversion functions are used to allow for all angles, including angles in the form of π/2 + n*π.

The following code demonstrates the two basic approaches, using the classic HP 32 (and which can be used on HP 32SII, DM32, HP 33S).

HP 32S Code

The variables:

X = (initial) position of the spring. When the spring is at it’s natural or neutral position, X = 0. (SI: unit m)

V = (initial) velocity of the spring. If V<0, the spring is getting shorter. If V>0, the spring is getting longer.

M = mass of the object attached to the end of the string (SI: unit kg)

K = spring’s force constant (N/m, kg/s^2)

A = amplitude of the spring (m)

Labels:

LBL I: Initialization, calculate W = √(K / M), and A = √(X^2 + V^2/W^2)

(HP 32S Memory: I: 34.5 bytes)

I01 LBL I

I02 INPUT K

I03 INPUT M

I04 INPUT X

I05 INPUT V

I06 RAD

I07 RCL K

I08 RCL÷ M

I09 SQRT

I10 STO W

I11 VIEW W

I12 RCL X

I13 x^2

I14 RCL V

I15 x^2

I16 RCL W

I17 x^2

I18 ÷

I19 +

I20 SQRT

I21 STO A

I22 VIEW A

I23 RTN

Sine Approach: x(t) = A * sin(W*T + B); W = √(K / M)

LBL S: Calculate B, then calculates X from time = 0 seconds to T seconds (per second)

(HP 32S Memory: S: 21.0 bytes, U: 16.5 bytes)

S01 LBL S

S02 RCL X

S03 RCL V

S04 RCL W

S05 ÷

S06 y,x→Θ,r

S07 x<>y

S08 STO B

S09 VIEW B

S10 INPUT T

S11 3

S12 10^x

S13 ÷

S14 STO T

U01 LBL U

U02 RCL T

U03 IP

U04 RCL× W

U05 RCL+ B

U06 SIN

U07 RCL× A

U08 STOP

U09 ISG T

U10 GTO U

U11 RTN

Cosine Approach: x(t) = A * cos(W*T + D); W = √(K / M)

LBL C: Calculate D

LBL D: Calculate X from time = 0 seconds to T seconds (per second)

(HP 32S Memory: C: 22.5 bytes, U: 16.5 bytes)

C01 LBL C

C02 RCL V

C03 RCL W

C04 ÷

C05 +/-

C06 RCL X

C07 y,x→Θ,r

C08 x<>y

C09 STO D

C10 VIEW D

C11 INPUT T

C12 3

C13 10^x

C14 ÷

C15 STO T

E01 LBL E

E02 RCL T

E03 IP

E04 RCL× W

E05 RCL+ D

E06 COS

E07 RCL× A

E08 STOP

E09 ISG T

E10 GTO E

E11 RTN

Examples

The calculator is set to FIX 5 mode.

Example 1:

K = 40, M = 1.33, V = 0, X = 0.3. Up to 5 seconds

Results:

W = 5.48408, A = 0.30000

	Sine Approach (LBL S)	Cosine Approach (LBL C)
Time (seconds)	B = 1.57080 (π/2)	D = 0
0	0.30000	0.30000
1	0.20921	0.20921
2	-0.00822	-0.00822
3	-0.22067	-0.22067
4	-0.29955	-0.29955
5	-0.19711	-0.19711

Example 2:

K = 48, M = 3, V = 0, X = -0.5. Up to 5 seconds

Results:

W = 4.00000, A = 0.50000

	Sine Approach (LBL S)	Cosine Approach (LBL C)
Time (seconds)	B = -1.57080 (-π/2)	D = 3.14159 (π)
0	-0.50000	-0.50000
1	0.32682	0.32682
2	0.07275	0.07275
3	-0.42193	-0.42193
4	0.47883	0.47883
5	-0.20404	-0.20404

Sources

“Simple Harmonic Motion” Hyperphysics. http://hyperphysics.phy-astr.gsu.edu/hbase/shm.html Retrieved February 7, 2025.

LibreTexts Physics. “15.2: Simple Harmonic Motion” https://phys.libretexts.org/Bookshelves/University_Physics/University_Physics_(OpenStax)/Book%3A_University_Physics_I_-_Mechanics_Sound_Oscillations_and_Waves_(OpenStax)/15%3A_Oscillations/15.02%3A_Simple_Harmonic_Motion Retrieved February 6, 2025.

Rodicek, Danny. Mathematics & Physics for Programmers Hingham, MA. 2005. ISBN 1-58450-330-0

Eddie

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

The content on this blog is 100% generated by humans. The author does not use AI engines and never will.

Casio Classpad fx-CP400: Solving a System of Linear Differential Equations Using Laplace Transforms

Casio Classpad fx-CP400: Solving a System of Linear Differential Equations Using Laplace Transforms

Introduction

The following procedure should work with the Casio Classpad family of calculators (300, 330, fx-CP400, fx-CP500). 

An attempt to use the Casio Classpad to solve the system of differential equations:

dxf/dt = A * xf(t) + B * yf(t) + g1(t)

dyf/dt = C * xf(t) + D * yf(t) + g2(t)

Initial conditions:  xf(0) = xf0,   yf(0) = yf0

The independent variable in all the functions is t.  The Casio Classpad's CAS functions laplace and invLaplace are used to accomplish this task.

Syntax

Laplace Transform of f(t):

laplace(f(t), t, s)

Laplace Transform of a Differential Equation:

laplace(diff eq, independent variable, dependent variable, s)

Inverse Laplace Transformation:

invLaplace(ℒ(s), s, t)

t:  parameter of the original function f(t)

s:  parameter of the transformed function ℒ(s)

The definition of a Laplace Transformation:

ℒ(f(t))(s) = ∫( f(t) * e^(-s*t) dt, t = 0 to ∞)

Now, even though we are working a system of differential equations, the laplace and invLapace commands can only work on function at a time.  Any other variable will be treated as a constant, even it was meant to represent another functions.  The substitution steps will account for this.  

The variables used to respect the Classpad's system variables.  

Here are the steps I used to solve:

dxf/dt = A * xf(t) + B * yf(t) + g1(t)

dyf/dt = C * xf(t) + D * yf(t) + g2(t)

Initial conditions:  xf(0) = xf0,   yf(0) = yf0

Spaces are added to readability.  

1.  Store dxf/dt in the variable f1.   Use xf' to represent dxf/dt.  

2.  Store dyf/dt in the variable f2.   Use yf' to represent dyf/dt.

3.  Execute the following and substitutions:

laplace( f1, t, xf, s ) ⇒ g1

g1 | xf(0) = xf0 and  yf / s = Mp ⇒ g1

4.  Execute the following and substitutions:

laplace( f2, t, yf, s ) ⇒ g2

g2 | yf(0) = yf0 and Lp = Mp ⇒ g2

g2 | xf /s = Lp ⇒ g2

Note Lp = ℒ( xf(t) ) and Mp = ℒ( yf(t) )

5.  Now solve the system of transformed equations: 

solve( {g1, g2}, {Lp, Mp} ) ⇒ lists

6.  Take the inverse laplace transforms:

invLaplace( getRight( lists[1], s, t ) ) ⇒ h1

invLaplace( getRight( lists[2], s, t ) ) ⇒ h2

The solutions are:

xf(t) = h1

yf(t) = h2

The getRight command extracts the right side of an equation.

Examples

Example 1:

xf ' = -0.08 * xf + 0.02 * yf + 6

yf ' = 0.08 * xf - 0.08 * yf

Initial conditions:  xf(0) = 0, yf(0) = 150

Example 2:

xf ' = 2 * xf + yf + t^2

yf ' = -2 * xf + 5

Initial conditions:  xf(0) = 2, yf(0) = 5

Casio Classpad Program:  laplaced

Download here:  https://drive.google.com/file/d/1zxTEynCjXmBrPU9pYqmXo3kF_so7fiI-/view?usp=share_link

Code:

' setup

SetRadian

' local variables

Local str1, str2, f1, f2

Local g1, g2, lists, h1, h2

Local Mp, Lp, t, s, xval, yval

' main

ClrText

InputStr str1, "xf'(xf,yf,g1(t))="

InputStr str2, "yf'(xf,yf,g2(t))="

Input xval, "xf(0)? "

Input yval, "yf(0)? "

StrJoin "xf'=", str1, str1

StrJoin "yf'=", str2, str2

strToExp(str1) ⇒ f1

strToExp(str2) ⇒ f2

' transform and solve

laplace(f1, t, xf, s) ⇒ g1

g1 | xf(0) = xval and yf/s = Mp ⇒ g1

laplace(f2, t, yf, s) ⇒ g2

g2 | yf(0) = yval and Lp = Mp ⇒ g2

g2 | xf/s = Lp ⇒ g2

Print "Laplace Transforms:"

PrintNatural g1

Print g1

PrintNatural g2

Print g2

solve({g1, g2},{Lp, Mp})⇒lists

invLaplace(getRight(lists[1]),s,t) ⇒ h1

invLaplace(getRight(lists[2],s,t) ⇒ h2

Print "Solutions: "

Print "xf = "

PrintNatural h1

Print h1

Print "yf = "

PrintNatural h2

Print h2

Note:  PrintNatural displays the expression in textbook form in a popup box, but it does not add anything to the text output terminal.  This is why I included both PrintNatural and Print commands.  

Source

Kreyszig, Erwin.  Advanced Engineering Mathematics  8th Edition  John Wiley & Sons, Inc:  New York, NY.  1999.  ISBN 0-471-15496-2

All original content copyright, © 2011-2023.  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 - Lambda Week: Solving Differential Equations with Runge Kutta 4th Order Method

Python - Lambda Week: Solving Differential Equations with Runge Kutta 4th Order Method

Welcome to Python Week!  This we we're going to cover calculus and the keyword lambda.

Note:  All Python scripts presented this week were created using a TI-NSpire CX II CAS.   As of June 2022, the lambda keyword is available on all calculators (in the United States) that have Python.   If you are not sure, please check your calculator manual. 

Solving Differential Equations

This following script solves the differential equation:

y' = dy/dx = f(x,y)

with initial condition y(x0) = y0

Repeat the steps for each step size h:

f1 = h * f(x0, y0)

f2 = h * f(x0 + h/2, y0 + f1/2)

f3 = h * f(x0 + h/2, y0 + f2/2)

f4 = h * f(x0 + h, y0 + f3)

x0 = x0 + h   (update x)

y0 = y0 + (f1 + 2*f2 + 2*f3 + f4)/6   (update y)

The small h is, the more accurate the calculated coordinates are.  

rk4lam.py:  Runge Kutta 4th Order Method

All answers are stored in the nested list t.  

from math import *

print("Runge Kutta 4th Order")

print("Math Module imported")

f=eval("lambda x,y:"+input("dy/dx = "))

# must call for float numbers one at a time

x0=eval(input("x0 = "))

y0=eval(input("y0 = "))

h=eval(input("h = "))

# ask for an integer

n=int(input("number of steps: "))

# set up table

t=[[x0,y0]]

# main loop

for i in range(n):

  f1=h*f(x0,y0)

  f2=h*f(x0+h/2,y0+f1/2)

  f3=h*f(x0+h/2,y0+f2/2)

  f4=h*f(x0+h,y0+f3)

  x0=x0+h

  y0=y0+(f1+2*f2+2*f3+f4)/6

  print([x0,y0])

  t.append([x0,y0])

print("Done.  Recall t for table.")

Examples

Results are rounded to five digits.  

Example 1:

dy/dx = 2*x*y + x,  y(0) = 0, h = 0.1, 5 steps

(Real solution:  y = 1/2 * (e^(x^2) - 1))

Results (which matches the exact results):

x = 0.1, y ≈ 0.00503

x = 0.2, y ≈ 0.02041

x = 0.3, y ≈ 0.04709

x = 0.4, y ≈ 0.08676

x = 0.5, y ≈ 0.14201

Example 2:

dy/dx = ln x + y, y(10) = 1

(Real Solution:  y = [∫(ln t * e^(-t) dt, t = 10 to x) + e^(-10)] * e^x

Exact Results:

x = 11, y ≈ 6.74551

x = 12, y ≈ 22.51732

x = 13, y ≈ 65.53659

x = 14, y ≈ 182.60824

x = 15, y ≈ 500.96552

Runge Kutta with h = 1, 5 steps:

x = 11, y ≈ 6.71066

x = 12, y ≈ 22.33376

x = 13, y ≈ 64.78988

x = 14, y ≈ 179.90761

x = 15, y ≈ 491.80768

Runge Kutta with h = 0.1, 50 steps:

x = 11, y ≈ 6.74551   (recall t[10])

x = 12, y ≈ 22.51728   (t[20])

x = 13, y ≈ 65.53643   (t[30])

x = 14, y ≈ 182.60766  (t[40])

x = 15, y ≈ 500.96358  (t[50])

This ends Python week for now, I hope you find this week helpful and resourceful.

Until next time,

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. 

Working with Laplace Transforms

Working with Laplace Transforms


Using Laplace Transforms to solve various differential equations


Table of Laplace Transforms

Function:  y(t)
Transformed Function:  f(s)

f(s) = ℒ( y(t) )
= ∫ y(t) ∙ exp(-s ∙ t) dt from t = 0 to t = ∞,  s is a complex number

a1, a2, and a are constants

ℒ( a1 ∙ y1(t) + a2 ∙ y2(t) ) = a1 ∙ f1(s) + a2 ∙f2(s)

ℒ ( y'(t) ) = s ∙ f(s) - y(0)

ℒ( t ∙ y(t) ) = -f'(s)

ℒ( 1 ) = 1 ÷ s

ℒ( a ) = a ÷ s

ℒ( t^n ) = n! ÷ (s^(n+1))

ℒ( t^n ÷ n! ) = 1 ÷ (s^(n+1))

ℒ( e^(a ∙ t) ) = 1 ÷ (s - a)

ℒ( sin(a ∙ t) ) = a ÷ (s^2 + a^2)

ℒ( cos(a ∙ t) ) = s ÷ (s^2 + a^2)


All the problems presented today will have the initial condition y(0) = 0.  



Problem 1:  y'(t) + a ∙ y(t) = b with y(0) = 0

Note: y'(t) = dy/dt.   All of the problems this weekend will have an initial condition.

 y'(t) + a ∙ y(t) = b

Applying the Laplace Transform:

s ∙ f(s) - y(0) + a ∙ f(s) = b ÷ s

f(s) ∙ (s + a) = b ÷ s

f(s) = b ÷ (s ∙ (s + a))



Partial fractions:

b ÷ (s ∙ (s + a)) = G ÷ s + H ÷ (s + a)

Simplifying and considering the numerator:

b + 0 ∙ s = a ∙ G + s ∙ (G + H)

Then:

b = a ∙ G
0 = G + H

G = b/a,  H = -G, H = -b/a


f(s) = (b/a) ÷ s - (b/a) ÷ (s + a)

Reverse the Laplace Transform:

y(t) = b/a - b/a ∙ e^(a ∙ t)

y(t) = b/a ∙ (1 - e^(a ∙ t))



Problem 2:  y'(t) + a ∙ y(t) = e^(b ∙ t) with y(0) = 0

y'(t) + a ∙ y(t) = e^(b ∙ t)

Applying the Laplace Transform:

s ∙ (f(s) - y(0) + a ∙ f(s) = 1 ÷ (s - b)

(s + a) ∙ f(s) = 1 ÷ (s - b)

f(s) = 1 ÷ ((s - b) ∙ (s + a))


Partial fractions:

1 ÷ ((s - b) ∙ (s + a)) = E ÷ (s - a) + F ÷ (s + a)

Simplifying and considering the numerator:

1 = E ∙ (s + a) + F ∙ (s - b)

0 ∙ t + 1 = E ∙ (s + a) + F ∙ (s - b)

Solving for E and F:

1 = a ∙ E - b ∙ F
0 = E + F

E = -F;  F = -1 ÷ (a + b), E = 1 ÷ (a + b)


f(s) = (1/(a + b)) ÷ (s - b) - (1/(a - b)) ÷ (s + a)

Reverse the Laplace Transform:

y(t) = e^(b ∙ t) ÷ (a + b) - e^(-a ∙ t) ÷ (a + b)

y(t) = (e^(b ∙ t) - e^(-a ∙ t)) ÷ (a + b)



Problem 3:  y'(t) + a ∙ y(t) = b ∙ t + c with y(0) = 0

y'(t) + a ∙ y(t) = b ∙ t + c

Applying the Laplace Transform:

s ∙ f(s) - y(0) + a ∙ f(s) = (b ∙ 1!) ÷ (s^(1 + 1)) + c ÷ s

s ∙ f(s) + a ∙ f(s) = b ÷ s^2 + c ÷ s

f(s) = (b + c ∙ s) ÷ (s^2 ∙ (s + a))


Partial fractions:

(b + c ∙ s) ÷ (s^2 ∙ (s + a)) = G ÷ (s + a) + (s ∙ H + J) ÷ s^2

= G ÷ (s + a) + H ÷ s + J ÷ s^2

Simplifying and considering the numerator:
 
b + c ∙ s + 0 ∙ s^2  = G ∙ s^2 + H ∙ s^2 + H ∙ a ∙ s + J ∙ s + J ∙ a

b = a ∙ J
c = a ∙ H + J
0 = G + H

H = -G; J = b/a, H = (a ∙ c - b) ÷ a^2, G = (b - a ∙ c) ÷ a^2


f(s) = ((b - a ∙ c) ÷ a^2) ÷ (s + a) + ((a ∙ c - b) ÷ a^2) ÷ s + (b/a) ÷ s^2

Reverse the Laplace Transform:

y(t) = ((b - a ∙ c) ÷ a^2) ∙ e^(-a ∙ t) + ((a ∙ c - b) ÷ a^2)  + (b/a) ∙ t





Problem 4:  y'(t) + a ∙ y(t) = b ∙ t^2 with y(0) = 0

y'(t) + a ∙ y(t) = b ∙ t^2

Applying the Laplace Transform:

s ∙ f(s) - y(0) + a ∙ f(s) = 2 ∙ b ÷ s^3

f(s) = (2 ∙ b) ÷ (s^3 ∙ (s + a))


Partial fractions:

(2 ∙ b) ÷ (s^3 ∙ (s + a)) = (F ∙ s^2 + G ∙ s + H) ÷ s^3 + K ÷ (s + a)

(2 ∙ b) ÷ (s^3 ∙ (s + a)) = F ÷ s + G ÷ s^2 + H ÷ s^3 + K ÷ (s + a)

Simplifying and considering the numerator:

0 ∙ s^3 + 0 ∙ s^2 + 0 ∙ s + 2 ∙ b = F ∙ s^3 + G ∙ s^2 + H ∙ s + F ∙ a ∙ s^2 + G ∙ a ∙ s + H ∙ a + K ∙ s^3

2 ∙ b = a ∙ H
0 = H + a ∙ G
0 = G + a ∙ F
0 = F + K

F = -K; K = (-2 ∙ b) ÷ a^3 ; F = (2 ∙ b) ÷ a^3; H = (2 ∙ b) ÷ a; G = (-2 ∙ b) ÷ a


f(s) = (2 ∙ b)/a^3 ÷ s - (2 ∙ b)/a^2 ÷ s^2 + (2 ∙ b)/a ÷ s^3 - (2 ∙ b)/a^3 ÷ (s + a)

I use Wolfram Alpha and CAS calculators to help verify whether the work is done correctly.  Sometimes I have deal with problems in parts.  

Reverse the Laplace Transform:

y(t) = (2 ∙ b)/a^3 - (2 ∙ b)/a^2 ∙ t + (b/a) ∙ t^2 - (2 ∙ b)/a^3 ∙ e^(-a ∙ t)

Note:  ℒ( t^n ÷ n! ) = 1 ÷ (s^(n+1)), ℒ( t^2 ÷ 2) = 1 ÷ (s^2)


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.