Solving initial value problems in Excel

Syntax

=IVSOLVE(eqns, vars, domain, [options])


Use IVSOLVE to solve an initial-value ordinary (ODE) or differential-algebraic (DAE) equation system stated in either of the canonical forms below: (the system can have as many equations as needed.)

An explicit ODE or DAE system

Example: dy1dt=f1(t,y1,y2,y3) dy2dt=f2(t,y1,y2,y3) With optional one or more algebraic equations: 0=g1(t,y1,y2,y3)

An implicit system coupled by a mass matrix

Example: m11m12m21m22y˙1(t) =f1(t,y1,y2)y˙2(t)=f2(t,y1,y2)


Initial Values

To complete description of the problem, initial values for the variables yi are required.

Required Inputs

eqns a reference to the system formulas (f1, f2,..,g1,..).

For a DAE system, algebraic equations must be listed after differential equations.

vars a reference to the system variables in the following order (t, y1, y2,.. ). Each state variable must be assigned an initial value.

For a DAE system, algebraic variables should be orderd last in vars.

domain a vector defining the start and end values for the time interval as well as the desired solution output time values.

FormatRemarks
{ti, tf}The interval [ti, tf] is divided uniformly according to the number of allocated rowsin the output array.
{{ti, tf, ndiv}The interval [ti, tf] is divided uniformly into ndiv subintervals.The allocated output array must have sufficient rows to hold the result.
{ti, t1, t2, .., tf}Solution results are reported for the supplied values only. The allocated output array must have sufficient rows to hold the result.

m the number of algebraic constraints for an explicit DAE system, or the full mass Matrix for an implicit system.

tol controls the absolute and relative error tolerances. Default values are 1.0e-4 for relative tolerance and 1.0e-6 for absolute tolerance. Various formats are permitted as described below.

FormatRemarks
rtolRelative tolerance. Absolute tolerance value is set to rtol/100. Fixed for all variables.
{rtol, atol}Relative and absolute tolerance values. Fixed for all variables
vector of rtolicustom relative tolerance for each system variable. Absolute tolerance is set to rtoli/100.
two-column array of {rtoli, atoli}custom relative and absolute tolerances for each system variable.

ctrl a set of key/value pairs for algorithmic control as detailed below.

Description of key/value pairs for algorithmic control

KeyALGOR
Admissible Values (String)RADAU5, BDF, ADAMS
Default ValueRADAU5
RemarksSee Algorithms Section for details.
KeyMAXSTEPS
Default Value100000 (Integer)
RemarksSets an upper bound on the maximum number of integration steps that can be carried out.

KeyINITSTEP
Default Value1.0e-5 (Real)
RemarksThe step size is quickly adapted. A small value is recommended if the system has fast initial transients.

KeySTEPLIMIT
Default ValueUnbounded (Real)
RemarksBy setting a limit, accuracy may be improved at the expense of speed.
KeyJACSTEP
Admissible Values (Real)0 < JACSTEP < 0.1
Default ValueThe default step is computed dynamically based on machine accuracy, preset limits and function metric. For an order(1) function it is approximately 1.0e-8.
RemarksThe default value generally produces accurate approximations; however relaxing the default value may aid convergence on some problems with unknown Jacobian such as parameterized problems.

KeyJACSCHEME
Admissible Values (Integer)

1 for First order Euler forward scheme

2 for Second order Euler forward scheme

Default Value1
RemarksFirst order scheme is generally sufficient.

jac a reference to the system analytical Jacobian matrix for the functions (f1, f2,..,g1,..) with respect to the variables (y1, y2,.. ).

For a system with two functions, the Jacobian matrix has the following order: f1y1f1y2f2y1f2y1

IVSOLVE is evaluated as an array formula in an allocated range of the spreadsheet. It computes and displays a formatted solution as illustrated in Figure 1 for a system with two state variables and a time domain [0 ,1]. The first column displays the time points, and the subsequent columns display the corresponding values of the state variables.

Figure 1: IVSOLVE formatted output
ABC
1ty1y2
20y1(0)y2(0)
30.1
40.2
..
121y1(1)y1(1)

When allocating a range for BVSOLVE, the following rules apply:

  • The number of allocated columns should match the number of variables.

  • The number of allocated rows is arbitrary. It determines the increment for the time interval subdivisions (only for display purpose).
    You can customize the output time points using any of the following formats in argument number 3.

    FormatRemarks
    {ti, tf}The interval [ti, tf] is divided uniformly according to the number of allocated rowsin the output array.
    {{ti, tf, ndiv}The interval [ti, tf] is divided uniformly into ndiv subintervals.The allocated output array must have sufficient rows to hold the result.
    {ti, t1, t2, .., tf}Solution results are reported for the supplied values only. The allocated output array must have sufficient rows to hold the result.

The following rate equations describe a chemical reaction

dy1dt=-0.04y1+104y2y3 dy2dt=0.04y1-104y2y3-3*107y22 dy3dt=3*107y22

This system time domain is [0,1000] with initial conditions y1 = 1, y2 = y3 = 0

Solution

We pick T1 to represent the time variable and Y1, Y2, Y3 to represent the differential (state) variables, and define the system right-hand-side formulas in terms of our selected variables in range A1:A3 as shown in Table 1. We also assign the initial conditions for the state variables as shown in Table 1.

Table 1
AY
1=-0.04*Y1+10000*Y2*Y31
2=0.04*Y1-10000*Y2*Y3-30000000*Y2^20
3=30000000*Y2^20

Next, we evaluate the array formula =IVSOLVE(A1:A3, (T1,Y1:Y3), {0,1000}) in range G4:J25 to obtain the solution shown in Table 2. The first argument to IVSOLVE is the reference to the system RHS formulas defined in A1:A3. The second parameter is the reference to the selected system variables. Note that we have used Excel union operator, (), to combine the time variable T1 and the state variables Y1:Y3 into one reference. The third parameter is the time domain for the problem which is passed using Excel constant array syntax, {}, for convenience. Since we have used default format for output, IVSOLVE reports the solution at a uniform time subdivision of the interval [0,1000] based on the number of allocated rows in the allocated range for output. You can control the output by allocating a larger range or using other optional formats for parameter 3 to control the subdivison or the exact outputed time points.

Table 2
GHIJ
4T1Y1Y2Y3
50100
6500.6932988.34608E-060.307238
71000.6176246.15491E-060.382913
81500.5706035.12544E-060.429936
92000.53624.48885E-060.464339
102500.5090244.04263E-060.491515
113000.4865733.70661E-060.513967
123500.4674643.44094E-060.533076
134000.4508313.22379E-060.549709
144500.4361313.04165E-060.56441
155000.4229752.88612E-060.577566
165500.4110622.75081E-060.589479
176000.4002082.63192E-060.600333
186500.3902442.52635E-060.610297
197000.3810272.43154E-060.619514
207500.3724682.34583E-060.628073
218000.3644912.26794E-060.63605
228500.3570232.19669E-060.643519
239000.352.13113E-060.650541
249500.3433822.07052E-060.65716
2510000.3371272.01432E-060.663415

Consider the following differential algebraic system of equations: dy1dt=-0.04y1+104y2y3 dy2dt=0.04y1-104y2y3-3*107y22 0=y1+y2+y3-1

This is a variation of Example 1 where y3 is treated as an algebraic variable. This system time domain is [0,1000] with initial conditions y1 = 1, y2 = y3 = 0

Solution

The solution steps are similar to Example 1. The RHS formulas are defined in A4:A6 as shown in Table 1 in terms of the selected variables T1 and Y1, Y2, Y3. Note that Y3 represents now an algebraic variable rather than a differential variable.

Table 1
A
4=-0.04*Y1+10000*Y2*Y3
5=0.04*Y1-10000*Y2*Y3-30000000*Y2^2
6=Y1+Y2+Y3-1

Next, we evaluate the array formula =IVSOLVE(A4:A6, (T1,Y1:Y3), {0,1000}, 1) in an allocated range to obtain the solution. Note that we pass 1 in the 4rth parameter to tell the IVSOLVE that we have one algebraic equation in the system. Algebraic equations and variables must be ordered last. Here the solver will treat the last formula as an algebraic equation along with variable Y3.

The new results, are only slightly different from Example 1, however, we extract below a sample to show that the solution satisfies the algebraic constraint Y1+Y2+Y3-1=0 as shown in Table 2.

Table 2
Y1Y2Y3Value of
Y1+Y2+Y3-1
1000
0.6928798.34419E-060.307113-3.1E-13
0.6172456.15373E-060.3827493.31E-13
0.5702325.12416E-060.429763-2.1E-13
0.5358464.48774E-060.4641497.77E-14
0.5086954.04186E-060.4913010

Next, we demonstrate how to use additional optional arguments.

Passing System Jacobian, and Custom Tolerances

We define the system Jacobian in M3:O5 and custom tolerances for variables Y1, Y2 and Y3 in the vector Q3:Q5 as shown in Table 3. We name these two ranges in Excel as jacob and rtol respectively.

Table 3
MNO Q
3-0.04=10000*Y3=10000*Y20.00001
40.04=-10000*Y3-30000000*2*Y2=-10000*Y20.0000001
51110

We evaluate the updated array formula =IVSOLVE(A4:A6, (T1,Y1:Y3), {0,1000,16}, 1, rtol, jacob) in an allocated array M7:P25 and obtain the new solution shown in Table 4. Notice that we requested in argument 3 to use exactly 16 subdivisions for the output time interval. This fills rows 8 through 24. Any extra rows in the allocated array are filled with zeros.

Table 4
MNOP
7T1Y1Y2Y3
80100
962.50.6692277.57291E-060.330766
101250.5915945.56659E-060.408401
11187.50.543624.62411E-060.456375
122500.5086854.04171E-060.491311
13312.50.4811973.63374E-060.518799
143750.458563.3264E-060.541437
15437.50.4393423.08361E-060.560655
165000.4226722.88522E-060.577325
17562.50.4079692.71896E-060.592028
186250.3948372.57688E-060.60516
19687.50.3829862.45358E-060.617012
207500.3722012.3452E-060.627796
21812.50.3623132.24888E-060.637685
228750.3531972.16257E-060.646801
23937.50.3447452.08461E-060.655253
2410000.3368752.01371E-060.663122
250000

The following system with 12 rate equations arises from chemical kinetics.

dy1dt=-0.1y1 dy2dt=0.1y1  +1500y4  +1500y5-50y2y3-100y2y12-10y2 dy3dt=10y2-0.1y3  -50y2 y3-50y10y3+1500y4 +1500y6 dy4dt=50y2y3-1502.5y4 dy5dt=100y2y12-1502.5y5 dy6dt=50y10y3-1502.5y6 dy7dt=100y10y12-1502.5y7 dy8dt=50y10-1505y8 dy9dt=2.5y4+2.5y5 +2.5y6+2.5y7 dy10dt=0.1y3+1500y6+1500y7+1500y8-50y10y3-100y10y12-10y10-50y10 dy11dt=5y8 dy12dt=10y10+1500y5+1500y7-100y2y12-100y10y12

We solve this system in the time interval [0,100] starting from initial conditions: y1 = 1 and 0 for the rest.

Solution

The steps to solve this system are similar to Example 1. We define the RHS system formulas in A26:A37 using variables T1, and Y1:Y12 as shown in Table 1. Variable Y1 is assigned the value 1 for the initial condition whereas the remaining variables Y2:Y12 are left blank consistent with the initial condition of 0.

Table 1
A
26=-0.1*Y1
27=0.1*Y1+1500*Y4+1500*Y5-50*Y2*Y3-100*Y2*Y12-10*Y2
28=10*Y2-0.1*Y3-50*Y2*Y3-50*Y10*Y3+1500*Y4+1500*Y6
29=50*Y2*Y3-1502.5*Y4
30=100*Y2*Y12-1502.5*Y5
31=50*Y10*Y3-1502.5*Y6
32=100*Y10*Y12-1502.5*Y7
33=50*Y10-1505*Y8
34=2.5*Y4+2.5*Y5+2.5*Y6+2.5*Y7
35=0.1*Y3+1500*Y6+1500*Y7+1500*Y8-50*Y10*Y3-100*Y10*Y12-10*Y10-50*Y10
36=5*Y8
37=10*Y10+1500*Y5+1500*Y7-100*Y2*Y12-100*Y10*Y12

To obtain the solution, we evaluate the array formula =IVSOLVE(A26:A37, (T1,Y1:Y12), {0,100}) in allocated array J4:V26. The results are plotted in the figure below.

IVSOLVE offers the following algorithms:

RADAU5
An implicit Runge-Kutta method of order 5 with adaptive step size and error control. Highly stable and suitable for nonlinear stiff problems. (Default algorithm.)
BDF
An implicit backward differentiation formula method (DLSODE) suitable for nonlinear stiff and non-stiff systems
ADAMS
An implicit Adams method (DLSODI) suitable for non-stiff systems.

  • Ernst Hairer and Gerhard Wanner, Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems (Springer Series in Computational Mathematics), 1996.
  • Alan C. Hindmarsh, ODEPACK, A Systematized Collection of ODE Solvers, in Scientific Computing, R. S. Stepleman et al. (Eds.), North-Holland, Amsterdam, 1983, pp. 55-64.
Question or Comment? Email us:
support @ excel-works.com
ExceLab: Transforming Excel into a Calculus Power House
ExceLab functions and methods are protected by USA Patents 9286286, 9892108, 10114812 and pending.
© 2015-2019, ExcelWorks LLC
Boston, USA