The Prothero-Robinson stiff ODE model

Bartel & Günther present a model for stiff systems as adapted from Kvaerno. The Prothero Robinson model that they use is defined by the following set of non-homogeneous equations:

 

The above equations represent sinusoidal solutions governed by the frequency w, with the stiffness factor governed by gamma. Bartel & Günther suggest -10E4 for the stiffness factor, with zero initial conditions at t=0, w=100 and eps =0.1. The solutions of the ODE’s are thus two variables with vastly differing settling rates, clearly justifying the classification as stiff. 

Plotted below is the solution obtained using the Generalized Rosenbrock Algorithm with the VELDS parameters.

 

The non-linear Coupled Pendulum

The equations for the coupled pendulum as presented in Hairer et al are: 

2^nd order ODE’s can be numerically integrated by converting into 1^st order equivalents. This can be performed by a straight-forward conversion, or using the State-Space approach. Once represented as 1^st order ODE’s, a numerical technique such as Runge-Kutta or Rosenbrock can be easily implemented to achieve the required output. Shown below is the solution to the coupled pendulum, integrated using the Generalized Rosenbrock method with the VELDS coefficients.

Veldhuizen’s Coefficients for Rosenbrock Methods (VELDS)

Veldhuizen developed two sets of coefficients for use with the Rosenbrock formulation. Shown below is one such set, “VELDS”, as presented in Hairer & Wanner’s ROS4 code.

              A21= 0.2000000000000000E01

       A31= 0.1750000000000000E01

       A32= 0.2500000000000000

       C21=-0.8000000000000000E01

       C31=-0.8000000000000000E01

       C32=-0.1000000000000000E01

       C41= 0.5000000000000000

       C42=-0.5000000000000000

       C43= 0.2000000000000000E01

       B1= 0.1333333333333333E01

       B2= 0.6666666666666667

       B3=-0.1333333333333333E01

       B4= 0.1333333333333333E01

       E1=-0.3333333333333333

       E2=-0.3333333333333333

       E3=-0.0000000000000000

       E4=-0.1333333333333333E01

       GAMMA= 0.5000000000000000

       C2= 0.1000000000000000E01

       C3= 0.5000000000000000

       D1= 0.5000000000000000

       D2=-0.1500000000000000E01

       D3=-0.7500000000000000

       D4= 0.2500000000000000

 

Used with the Generalized Rosenbrock formulation, these coefficients can be implemented fairly easily for numerical solution of 1^st order ODE’s.

The Lorenz Attractor

The Lorenz attractor  is a series of first order ODE’s that was first studied by Edward Lorenz. The ODE’s as adapted from Hairer et al, along with the initial conditions, were solved using the Rosenbrock technique.

 

 

 

Plotted below are the solutions to the equations. 

The Van der Pol Oscillator

Presented here is a numerical solution to the Van der PolOscillator. The form of the ODE’s has been chosen from Hairer et al. The equations are:

The Initial conditions y1(0)=-0.5, y2(0)=0, A = 2.008619861. The non-linearity in the equations is fairly clear from the plots below. The solution of y1 remains stable for several seconds between sudden changes in response, clearly following the behavior of y2.