Numerics Lab

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 u_{1}, u_{2}, u_{3}, and u_{4} as unknowns, the differential equations are du_{1}/dt = u_{3}, du_{2}/dt = u_{4}, du_{3}/dt = u_{1} + 2·u_{4} + b·(u_{1} + a)/D_{1} - a·(u_{1} - b)/D_{2}, and du_{4}/dt = u_{2} - 2·u_{3} - b·u_{2}/D_{1} - a·u_{2}/D_{2}, where D_{1} and D_{2} are both functions of u_{1} and u_{2}.

**Autocatalytic Scheme**

Solves a reactor dynamics problem involving an autocatalytic reaction scheme. With y_{1}, y_{2}, and y_{3} as unknowns, the differential equations are dy_{1}/dt = a - y_{1} - b·y_{1}·y_{3}, dy_{2}/dt = y_{1} - y_{2}·y_{3}, and y_{3}/dt + y_{2}·y_{3} - b·y_{1}·y_{3}. 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 y_{1} and y_{2} as unknowns, it solves the differential equations dy_{1}/dt = alpha - y_{1} - 4·y_{1}·y_{2} / (1 + y_{1}·y_{1}) and dy_{2}/dt = beta * y_{1} * (1 - y_{2} / (1 + y_{1}·y_{1})). 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 = (x^{2} + y^{2}) (1 + 2x + 5x^{2} + 6x^{3} + 6x^{4} + 4x^{5} + x^{6} - 3y^{2} - 2xy^{2} + 8x^{2}y^{2} + 8x^{3}y^{2} + 3x^{4}y^{2} + 2y^{4} + 4xy^{4} + 3x^{2}y^{4} + y^{6})

**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 C_{a}, C_{b}, and C_{c} as unknowns, the equations are dC_{a}/dt = -k_{1}·C_{a}, dC_{b}/dt = k_{1}·C_{a} - k_{2}·C_{b}, and dC_{c}/dt = k_{2}·C_{b}.

**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 = sin^{2}(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 q_{1}, q_{2}, v_{1}, v_{2}, and z as unknowns, the differential equations are dq_{1}/dt = v_{1}, dq_{2}/dt = v_{2}, dv_{1}/dt = -z·q_{1}, and dv_{2}/dt = -z·q_{2} - g, subject to the algebraic constraint 1 = q_{1}·q_{1} + q_{2}·q_{2}.

**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 q_{1}, q_{2}, v_{1}, v_{2}, and z as unknowns, the differential equations are dq_{1}/dt = v_{1}, dq_{2}/dt = v_{2}, dv_{1}/dt = -z·q_{1}/r, dv_{2}/dt = -z·q_{2}/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.