Solving boundary value problems in Excel

Syntax

=BVSOLVE(eqns, vars, bcs, domain, [options])


Use BVSOLVE to solve a boundary-value ODE or DAE system in the form: (the system can have as many equations as needed)

Example: dy1dx=f1(x,y1,y2,y3) dy2dx=f2(x,y1,y2,y3) With optional one or more algebraic equations: 0=g1(x,y1,y2,y3)

Boundary Conditions

To complete description of the problem, boundary conditions implicating the differential variables are required. General nonlinear boundary conditions can be specified at any location within the the system spatial domain [xs, xe].

Required Inputs

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

vars a reference to the system variables in the following order (x, y1, y2,.. ).

Optional. You can define starting guess for (y1, y2,.. ). The guess can be a constant or a function of x.

bcs A reference to a 2-column array defining the boundary conditions. The first column specifies boundary conditions points, and the 2nd column specifies corresponding boundary condition formulas. The boundary conditions points must be ordered from left to right.

Boundary condition points may be repeated if multiple boundary conditions are specified at the same location.

A boundary condition formula is arranged with respect to zero on one side. For Example, suppose that the boundary condition (BC) at x0 is y1= y2+1, then the BC formula may be defined as =Y1-Y2-1.

domain a vector defining the start and end values for the spatial domain. Additional optional data controls the output spatial values in the computed solution as described below.

FormatRemarks
{xs, xe}The domain [xs, xe] is divided uniformly according to the number of allocated rowsin the output array.
{{xs, xe, ndiv}The domain [xs, xe] is divided uniformly into ndiv subintervals.The allocated output array must have sufficient rows to hold the result.
{xs, x1, x2, .., xe}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 equations for a DAE system.

For a DAE system, algebraic equations and variables should be listed last in eqns and vars.

rtol relative error tolerance. you may specify a fixed value for all variables, or supply vector of tolerances for each variable in vars. Default is 1.0e-4.

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

Description of key/value pairs for algorithmic control

KeyMAXMESH
Default Value1000 (Integer)
RemarksSets an upper bound on the mesh size for the spatial domain.

KeyINITMESH
Default Value7 (Integer)
RemarksSets the initial mesh size which is refined adaptively during iterations.

KeyKORDER
Admissible Values (Integer)2 ≤ KORDER ≤ 7
Default Value4
RemarksNo. of collocation points per mesh subinterval. A higher value will increase computational cost.It may help with a highly oscillating system.

KeyLINEAR
Admissible Values (Boolean)True or False
Default ValueFalse
RemarksSet to TRUE if the problem is linear with linear boundary conditions

KeyERROREST
Admissible Values (Boolean)True or False
Default ValueFalse
Remarks

Set to true to report global error estimates for the differential variables in the solution output array.The last row of the allocated output array will be reserved for the error estimates.

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.

sysjac 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

bcjac a reference to the boundary conditions analytical Jacobian matrix with respect to the differential variables only (y1, y2,.. ).

BVSOLVE 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 spatial domain [0 ,1]. The first column displays the spatial points, and the subsequent columns display the corresponding values of the state variables.

Figure 1: BVSOLVE formatted output
ABC
1xy1y2
20y1(0)y2(0)
30.1
40.2
..
121y1(1)y1(1)
13Rel. Errorerror y1error y2

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 domain subdivisions (only for display purpose).
    You can customize the output spatial points using any of the following formats in argument number 5.

    FormatRemarks
    {xs, xe}The domain [xs, xe] is divided uniformly according to the number of allocated rowsin the output array.
    {{xs, xe, ndiv}The domain [xs, xe] is divided uniformly into ndiv subintervals.The allocated output array must have sufficient rows to hold the result.
    {xs, x1, x2, .., xe}Solution results are reported for the supplied values only. The allocated output array must have sufficient rows to hold the result.
  • If the control key ERROREST value is set to True (see examples below), the last row in the allocated output array will display the estimated maximum global error for each variable.

Consider the following 2nd-order nonlinear differential equation in the spatial domain [0,1]

ϵ u''=-euu'+12 πsin12 πxe2u,     0x1 with end points boundary conditions: u0=0,  u1=0

Solution

Using standard substitution technique we convert the 2nd-order equation to two 1st-order equations

dudx=z dzdx=1ε(-euz+12 πsin12 πxe2u)

Working with the named variables x, y, z, and eps for cells B2:B5, we define the system equations and boundary conditions as shown in Table 1 below. The colored ranges represent the required input arguments to BVSOLVE. They are: the RHS-side formulas in B7:B8, the selected system variables B2:B4 (named x, y, z,), and the boundary conditions matrix C2:D3.

Table 1
ABCD
1VariablesBCs
2x0=u
3u1=u
4z
5eps0.01
6Equations
7du/dx=z
8dz/dx=(-EXP(u)*z+0.5*PI()*SIN(0.5*PI()*x)*EXP(2*u))/eps

Next, we evaluate the array formula =BVSOLVE(B7:B8, (x,y,z), C2:D3, {0,1}, , , {"ERROREST", True}) in an allocated array G4:I25 and obtain the solution shown in Table 2 below.

  • We have used Excel union operator, (), to combine the system variables x,y,z into one reference for argument 2. Alternatively, we could pass the raw address B2:B4 instead.
  • We have used Excel constant array syntax, {}, to pass the spatial domain [0,1] in argument 4.
  • Note that we skip over optional parameters 5 and 6 and pass the key ERROREST in parameter 7 to report global error estimates in the last row of the output array.
The obtained result is plotted in Figure 1.

Table 2
GHI
4xuz
500-9.22395
60.05-0.3043492-3.93625
70.1-0.44424441-1.92551
80.15-0.51353621-0.95066
90.2-0.54647914-0.41601
100.25-0.55866051-0.09636
110.3-0.557967560.110589
120.35-0.54860750.256436
130.4-0.532874280.368873
140.45-0.512016770.463304
150.5-0.486695810.548502
160.55-0.457234490.629561
170.6-0.423757590.709506
180.65-0.386270750.790212
190.7-0.344702970.872914
200.75-0.298930930.958517
210.8-0.248790481.047768
220.85-0.194081621.14137
230.9-0.134569021.240042
240.95-0.069979941.344568
25101.455833
26Rel. Errors1.87103E-068.35E-05

Figure 1

Optional Parameters

For demonstration purpose, we show how to define and pass the optional system and boundary conditions jacobians, well as custom tolerances. (In practice, passing in analytical Jacobians is unnecessary as the internal numerical approximations by the solver is usually accurate and sufficient.) We define the system jacobian in D8:E9, and the boundary conditions jacobian in D11:E12 as shown in Table 3 below.
The jacobian definition must respect the order of state variables, u, z as specified in argument 2 for the solver.

Table 3
DE
7System Jacobian
801
9=-(EXP(u)*z-PI()*SIN(PI()*x/2)*EXP(2*u))/eps=-EXP(u)/eps
10BCs Jacobian
1110
1210

We evaluate the updated formula
=BVSOLVE(B7:B8, (x,y,z), C2:D3, {0,1}, 0, {0.0001,0.00001}, {"ERROREST",True}, D8:E9, D11:E12)
to obtain the new solution. Given that we have only 2 state variables, it is convenient to use an array constant to pass the tolerances vector in argument 6. The new results are essentially the same as the default results shown above in Table 2.

Consider the following stiff nonlinear 2nd-order differential equation in the domain [0,1]
ε Axyy''-2.42-εA'xyy'+y'y+A'xAx1-0.42y2=0,  0x1 Ax=1+x2 with end points boundary conditions 0=0.9129,        y1=0.375

Using standard substitution technique we convert the 2nd order equation to two 1st order equations

dydx=z dzdx=1ε Axy 2.42-εA'xyz-zy-A'xAx1-0.42y2

Solution

We solve the above stiff system for ε=0.01.

The solution steps are similar to Example 1. We work with the named variables, x, y, z, A, dA, eps corresponding to cells B2:B7 to define the system formulas and boundary condition matrix as shown in Table 1 below.

Table 1
ABCD
1VariablesBCs
2x0=y-0.9129
3y1=y-0.375
4z
5A=1+x^2
6dA=2*x
7eps0.01
8Equations
9dy/dx=z
10dz/dx=((2.4/2-eps*dA)*y*z-z/y-dA/A*(1-0.4/2*y^2))/(eps*A*y)

Next, we evaluate the array formula =BVSOLVE(B9:B10, (x,y,z), C2:D3, {0,1}, , , {"ERROREST", True}) in an allocated array H4:J26 to obtain the solution. However, the solver displays a fatal error message instead:

Explanation

An error message can often be quite helpful. The caption of the window indicates the cell and solver that has caused the error, and the message often conveys details about the specific argument responsible for the error.
In this case, the error is in argument 1 which holds the system equations. Clearly, the solver encountered division by zero while evaluating one of the system equations. A quick inspection shows that formula B10 involves division by the variable y which has a default starting guess of zero. To fix this error, we simply define a non-zero starting guess of 1 for y in cell B3. The solver now converges successfully and produces the following result (Table 2) which is plotted in Figure 1.

Table 2
HIJ
4xyz
500.91290.835223
60.050.9546460.834019
70.10.9962370.828936
80.151.0374680.819592
90.21.0781290.806201
100.251.1180260.789141
110.31.156990.7689
120.351.1948720.745969
130.41.2315420.72026
140.451.266760.684865
150.51.2986490.549927
160.551.305552-0.73972
170.61.083584-10.5522
180.650.603002-2.53229
190.70.547885-0.76692
200.750.511759-0.68353
210.80.479252-0.6188
220.850.449708-0.56442
230.90.422687-0.51745
240.950.397869-0.47612
2510.375-0.43931
26Rel. Errors8.33E-070.000237

Figure 1

This equation models a uniformly loaded beam of variable stiffness simply supported at both ends.

z''''=1-6x2z'''-6xz''x3 ,   1x2 z1=0,   z''1=0,   z2=0,   z''2=0

Using standard substitution technique we convert the 4th order equation to four 1st order equations.

Let u=z',         v=u'=z'',        w=v'=z''

dwdx=1-6x2w-6xvx3 dvdx=w dudx=v dzdx=u

The equivalent boundary conditions are z1=0,  v1=0,  z2=0,  v2=0

Solution

The solution procedure is similar to Example 1 above. Working with the named variables x,u,v,w and z for cells B2:B6, we define the system equations and boundary conditions as shown in Table 1 below.

Table 1
ABCD
1VariablesBCs
2x1=z
3u1=v
4v2=z
5w2=v
6z
7Equations
8dw/dx=(1-6*x^2*w-6*x*v)/x^3
9dv/dx=w
10du/dx=v
11dz/dx=u

To obtain the solution, We evaluate the array formula =BVSOLVE(B8:B11, B2:B6, C2:D5, {1,2}, , , {"ERROREST", True}) in allocated array H4:L26. The results are shown in Table 2 and plotted in Figure 1 below.

Table 2
HIJKL
4xwvuz
51-0.51.55242E-300.0171320
61.05-0.33011-0.020516090.0165840.000847162
71.1-0.20832-0.033809190.01520.001644528
81.15-0.12078-0.041916650.0132890.002358447
91.2-0.05787-0.046296290.0110710.002968346
101.25-0.0128-0.0480.0087040.003463059
111.30.019257-0.047792450.0063020.003838168
121.350.041773-0.046232790.0039470.004094076
131.40.057268-0.043731780.0016950.004234597
141.450.067583-0.04059207-0.000420.004265921
151.50.074074-0.03703704-0.002360.004195851
161.550.077746-0.03323151-0.004120.00403324
171.60.079346-0.02929688-0.005680.003787576
181.650.079432-0.02532209-0.007040.003468679
191.70.078423-0.02137187-0.008210.003086471
201.750.076635-0.01749271-0.009180.002650819
211.80.074303-0.01371742-0.009960.002171412
221.850.071605-0.01006851-0.010560.001657686
231.90.068677-0.00656072-0.010970.00111876
241.950.065617-0.00320302-0.011210.0005634
2520.06251.73472E-18-0.01129-1.0842E-19
26Rel. Errors3.17E-074.16988E-085.59E-098.15553E-10

Figure 1

Interior boundary conditions

We can easily study the effect of the beam support locations. We solve the same problem but with the beam supported at left end and center (x=1 and x=1.5) with the right end free. The only required change to the problem setup is to change the values of the boundary points in cells C4 and C5 in Tabe 1 from 2 to 1.5. We re-run the solver and obtain the new solution plotted below.

Consider the following differential algebraic (DAE) system of equations du1dx=ε+u2-sin(x)y+cos(x) du2dx=cos(x) du3dx=y 0=u1-sin(x)*(y-ex)

With boundary conditions

u10=sin0,   u21=sin1,   u30=1,   0x1

For small values of ε, this problem has the following solution

u1=u2= sinx;  u3=1; y=0 (See Ascher and Spiteri. Collocation software for boundary value differential-algebraic equations. SIAM Journal on Scientific Computing. Volume 15 938-952. 1995)

Solution

Working with the named variables x, u1_, u2_, u3_, y and eps for cells B2:B7, we define the system equations and boundary conditions as shown in Table 1 below.

Table 1
ABCD
1VariablesBCs
2x0=u1_-SIN(0)
3u1_0=u3_-1
4u2_1=u2_-SIN(1)
5u3_
6y
7eps0.001
8Equations
9du1/dx=(eps+u2_-SIN(x))*y+COS(x)
10du2/dx=COS(x)
11du3/dx=y
120=(u1_-SIN(x))*(y-EXP(x))

To obtain the solution, We evaluate the array formula =BVSOLVE(B9:B12, B2:B6, C2:D4, {0,1}, 1) in allocated array H4:L26 and obtain the results shown in Table 2 below.

  • The first argument B9:B12 is a reference to the system RHS formulas. Note that the algebraic equation B12 is ordered last as required for a DAE system.

  • The 2nd argument B2:B6 is a reference to the system variables. Alternatively, we could pass the equivalent defined names by combining them with Excel union operator as follows (x, u1_, u2_, u3_, y). Note that the algebraic variable B6 (y) is listed last as required for a DAE system.

  • We used Excel array constant syntax to pass the spatial domain [0, 1] in argument 4.

  • We used optional argument 5 to tell BVSOLVE that the last equation in B9:B12 is an algebraic equation.

Table 2
HIJKL
4xu1_u2_u3_y
500011.238E-05
60.050.0499790.04997916.752E-06
70.10.0998330.0998331-4.951E-07
80.150.1494380.1494381-1.035E-06
90.20.1986690.1986691-1.459E-06
100.250.2474040.24740414.503E-06
110.30.295520.295521-5.653E-06
120.350.3428980.3428981-1.840E-05
130.40.3894180.38941817.336E-06
140.450.4349660.4349661-4.392E-06
150.50.4794260.4794261-4.436E-05
160.550.5226870.52268715.825E-06
170.60.5646420.5646421-4.189E-07
180.650.6051860.6051861-8.164E-07
190.70.6442180.6442181-1.141E-06
200.750.6816390.68163913.401E-06
210.80.7173560.7173561-4.072E-06
220.850.751280.751281-1.321E-05
230.90.7833270.78332714.991E-06
240.950.8134160.8134161-2.817E-06
2510.8414710.8414711-2.832E-05

BVSOLVE is based on the COLDAE collocation method with adaptive mesh refinement.

  • U. Ascher and R. Spiteri. Collocation software for boundary value differential-algebraic equations, SIAM Journal on Scientific Computing. 15,938-952, (1994).
  • Ascher U M, Mattheij R M M and Russell R D. Numerical Solution of Boundary Value Problems for Ordinary Differential Equations, Prentice Hall, (1988).
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