Numerics Lab

Examples

Examples are provided as scripts that you can run by choosing the Example Scripts... item from the File menu.  Some of the current examples are described briefly below.

Astronomical Orbit

Solves the equations of motion in a plane for a small body under the influence of two large bodies of masses a and b = 1 - a. With u1, u2, u3, and u4 as unknowns, the differential equations are du1/dt = u3, du2/dt = u4, du3/dt = u1 + 2·u4 + b·(u1 + a)/D1 - a·(u1 - b)/D2, and du4/dt = u2 - 2·u3 - b·u2/D1 - a·u2/D2, where D1 and D2 are both functions of u1 and u2.

Autocatalytic Scheme

Solves a reactor dynamics problem involving an autocatalytic reaction scheme. With y1, y2, and y3 as unknowns, the differential equations are dy1/dt = a - y1 - b·y1·y3, dy2/dt = y1 - y2·y3, and y3/dt + y2·y3 - b·y1·y3.  When a = 3, the problem has a Hopf bifurcation at b = 1.30176.

Double Pendulum

One pendulum attached to another pendulum. A famous system illustrating chaotic behavior, formulated here as a high index problem.

Drawing Circles

Illustrates the dangers of applying numerical methods to the solution of differential equations. Using x and y as unknowns, it solves the equations dx/dt = -y and dy/dt = x.  With a constant step size of 0.02, the Implicit Midpoint method generates a circle, as expected, while the Backward Euler method generates a spiral.

Electric Circuit

A small circuit containing resistors, transistors, a capacitor, and a power source. It is described by a set of four algebraic equations and one differential equation.

Exponential Decay

Simple example with only one unknown, y. It solves the differential equation dy/dt = -y.

Hessenberg Three

Index 3 problem in Hessenberg form for x, y and z as the unknowns. It solves the differential equations dx/dt = y and dy/dt = z subject to the algebraic constraint x = 1 - exp(-2·t).

Hessenberg Two

Index 2 problem in Hessenberg form for x and y as the unknowns. It solves the differential equation dx/dt = -y subject to the algebraic constraint x = sin(2·pi·t).

Hopf Bifurcation

With y1 and y2 as unknowns, it solves the differential equations dy1/dt = alpha - y1 - 4·y1·y2 / (1 + y1·y1) and dy2/dt = beta * y1 * (1 - y2 / (1 + y1·y1)).  The problem is solved twice, for values of beta lower and higher than the critical value beta = 3·alpha/5 - 25/alpha.

Kepler Equation

Solves Kepler's Equation, M = E - e·sin(E) for an orbit with eccentricity e = 0.7 and period T = 365.

Long Integration

Illustrates the dangers of applying numerical methods to the solution of differential equations. Using y as unknown, it solves the equation dy/dt = y·cos(t), whose exact solution is y = exp(sin(t)). The solution obtained using the Forward Euler method with a constant step size of 0.1 deviates significantly from this exact solution.

Lorenz Attractor

Solves the famous problem illustrating chaotic behavior. As expected for chaotic systems, this problem is very sensitive to the choices you make of error tolerance, step size, and solution method.

Pear Curve

Solves the single nonlinear equation: 2 = (x2 + y2) (1 + 2x + 5x2 + 6x3 + 6x4 + 4x5 + x6 - 3y2 - 2xy2 + 8x2y2 + 8x3y2 + 3x4y2 + 2y4 + 4xy4 + 3x2y4 + y6)

Reactions in Series

Solves the mass balances in a batch reactor where a first order reaction turns A into B and another first order reaction turns B into C. With Ca, Cb, and Cc as unknowns, the equations are dCa/dt = -k1·Ca, dCb/dt = k1·Ca - k2·Cb, and dCc/dt = k2·Cb.

Robotic System

A high index problem. The robotic system is made up of two arms connected to each other, the first one fixed at the origin, while the second one is subject to the vertical constraint: y = sin2(wt). Fast changes in constraint forces make this a very difficult problem to solve.

Simple Pendulum

One of the best known index 3 problems, solves the equations of motion for a mass attached to a rigid rod. With q1, q2, v1, v2, and z as unknowns, the differential equations are dq1/dt = v1, dq2/dt = v2, dv1/dt = -z·q1, and dv2/dt = -z·q2 - g, subject to the algebraic constraint 1 = q1·q1 + q2·q2.

Spring Pendulum

When the rigid rod in the simple pendulum is replaced with a spring, the index of the problem is lowered to 1. With q1, q2, v1, v2, and z as unknowns, the differential equations are dq1/dt = v1, dq2/dt = v2, dv1/dt = -z·q1/r, dv2/dt = -z·q2/r - g, subject to the algebraic constraint e·z = r - 1, where r is the length of the spring and 1/e is the spring constant.

Stiff Problem

With y as the unknown, solves the equation dy/dt = -100·(y - sin(t)), whose solution varies rapidly initially but slows afterwards.