Differential Equations: Half-Increment Solution, Numerical Methods
11-27-2016, 06:41 PM
Post: #1
 Eddie W. Shore Senior Member Posts: 1,383 Joined: Dec 2013
Differential Equations: 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:

Free-Fall
d^2y/dt^2 = g
Set C as “9.80665” (SI) or “32.1740468” (US)

Free-Fall with Friction
d^2y/dt^2 = g - α (dy/dt)^2
(α = F/m)
Set C as “g - α * A^2”
(sub numeric values for g, α)

Spring
d^2x/dt = -k/m * x
Set C as “-k/m * T”
(sub numeric values for k, m)

Pendulum
d^2θ/dt = -α*sin(θ)
(α = -g/l)
Set C as “-α * sin(Y)”
(sub numeric values for α)

Damped, Driven Oscillations
d^2x/dt = -α*x – β*dx/dt + γ * sin(ω*t)
Set C as “-α*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.

Code:
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;

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