Solving Partial Differential Equations in Excel



Syntax

=PDASOLVE(eqns, vars, bcs, xdom, tdom, [options])


PDASOLVE is a versatile solver for partial differential equations that supports advanced modeling capabilities including:

  • Equations defined over multiple regions with discontinuous properties.
  • Flux continuity at interfaces between regions of discontinuous properties.
  • Flexible boundary conditions locations.
  • Flexible spatial domains definitions.
  • Algebraic constraints on variables.

Use PDASOLVE to solve a system of partial differential equations the following forms: (the system can have as many equations as needed)

An explicit PDE or PDA system

u1t=f1t,x,u,ux,uxx u2t=f2t,x,u,ux,uxx with optional algebraic equations 0=g1t,x,u,ux,uxx

An implicit system coupled by a mass matrix

m11m12m21m22u1/t=f1t,x,u,ux,uxxu2/t=f2t,x,u,ux,uxx

where u=[u1,u2], ux=[u1,x,u2,x], uxx=[u1,xx,u2,xx]

Initial Conditions

Initial value for each variable ui is required.

Boundary Conditions

PDASOLVE supports general boundary conditions that can be specified at arbitrary locations within the system's spatial domain. The boundary condition types are described below:

Dirichlet
A constant or time dependent condition imposed on a single differential variable.
Neumann
A constant or time dependent condition imposed on a single differential variable first derivative.
Robin
A general condition imposed on one or more differential variables and their first derivatives.
Continuity
A special condition to enforce continuity flux across an interface separating two regions with discontinuous material properties. A flux is defined as the product of a property such as conductivity and the first derivative of a differential variable.

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, x, u1, u2,.. , u1,x, u2,x,.., u1,xx, u2,xx,.. ). The state variables (u1, u2,..) must be assigned initial conditions. An initial condition can be a constant or function of the spatial variable x.

All the variables and derivatives must be defined in the order listed above in vars regardless if they are used in the equations or not.

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

bcs A reference to a 3-column matrix defining the system's boundary conditions. The first column specifies the x locations of the boundary conditions, the 2nd column specifies the types which can by any of the letters 'D', 'N', 'R' or 'C', and the 3rd column specifies the boundary condition formulas.

A boundary condition formula is arranged with respect to zero on one side.

D: Dirichlet
The formula must reference a single differential variable. For example, =u-1.
N: Neumann
The formula must reference a single differential variable first derivative. For example, =ux-1.
R: Robin
You can enter any arbitrary formula involving variables and/or their first derivatives.
Example, =k*ux-1, =ux-h*(u-tinf), =u1-u2^2*u1x.
C: Continuity
The formula should be an expression of the flux you want to enforce its continuity across the interface. For example, =k*ux.

Example of a boundary conditions matrix

EFG
100D=U1-1
111R=UX1-0.5*U1-2
120.5C=K1*UX1

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

If an equation is defined only over a sub region of the spatial domain, also use argument 7 to define custom regions.

FormatRemarks
{xs, xe}The domain [xs, xe] is divided uniformly according to the allocated output array size.
{xs, xe, ndiv}The domain [xs, xe] is divided uniformly into ndiv subintervals.
{xs, x1, x2, .., xe}Solution results are reported for the supplied values only.

tdom 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 allocated output array size.
{ti, tf, ndiv}The interval [ti, tf] is divided uniformly intondiv subintervals.
{ti, t1, t2, .., tf}Solution results are reported for the supplied values only.

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

rgns A 2-column matrix defining the subdomain for each equation. Each row specifies the left and right region endpoints for the corresponding equation in eqns.

. By default, all equations in the system are assumed to share the entire spatial domain defined in argument 4.

tol controls the absolute and relative error tolerances for temporal integration algorithm. 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

KeyNDRVOUT
Admissible Values (Integer)
  1. report only (u1, u2,..) in output.
  2. include first derivatives (u1,x, u2,x,..) in output.
  3. include first and second derivatives (u1,x, u2,x,..),(u1,xx, u2,xx,..) in output.
Default Value0
RemarksThe allocated results array must have sufficient columns to hold output.

KeyFORMAT
Admissible Values (String)

XCOL1
Column 1 of the results array will contain the output spatial points.A block of columns for the solution variables at each output temporal point will be reported in order following column 1.

TCOL1
Column 1 of the result array will contain the output temporal points.A block of columns for the solution variables at each output spatial point will be reported in order following column 1.

Default ValueXCOL1
RemarksXCOL1 format is convenient to plot snapshots of solution variables at a desired time values.
TCOL1 format is convenient to plot transient views of solution variables at a desired spatial locations.
See Output Format section for more detail.
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.
KeyMAX_GRID_SPACING
Default Value0.01 (Real)
RemarksSets a a limit on the maximum distance between adjacent grid nodes. All distances between adjacent grid nodes will be less than or equal toMAX_GRID_SPACING.
This value has direct impact on computational cost.

KeyMIN_GRID_SPACING
Default Value0.000001 (Real)
RemarksSets a limit on the minimum distance between adjacent grid nodes. All distances between adjacent grid nodes will greater than or equal toMIN_GRID_SPACING.

KeyMAX_GRID_NODES
Default Value2000 (Integer)
RemarksSets a limit on the maximum allowed number of generated nodes in the system spatial domain.

KeyMIN_GRID_NODES
Default Value10 (Integer)
RemarksSets a limit on the minimum number of generated nodes in any subregion of the system spatial domain.
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.

To illustrate PDASOLVE output layout, we consider a 2-equation system with the following variables (t, x, u1, u2, u1,x, u2,x, u1,xx, u2,xx)

Snapshot View Format

The snapshot view is the default format. It allows you to easily plot snapshot views for the variables at desired time points. In this format, the solution is laid out as shown in Table 1. Output spatial points are listed in the first column starting at row 3, whereas output time points are listed in blocks in the first row starting at column 2 with variables names listed in the second row.

Table 1
ABCDEFGH..
1tt0t0t1t1t2t2..tf
2xu1u2u1u2u1u2..u2
3xs
4x1u2(x1,t0)
5x2
....
Nxe
  • The spatial and time domain are divided uniformly according to the number of allocated rows and columns for the output array. You can control x and t values by specifying your own step or custom values in parameters 4 and 5 for PDASOLVE.

  • If you allocate a larger array than needed, unused rows and columns will be filled with zeros. You can find out the minimum array size needed by evaluating PDASOLVE formula initially in one cell.

We can request to report 1st and 2nd derivative variables using the key NDRVOUT with value 1 or 2. Table 2 shows the modified layout with NDRVOUT value set to 1, in which values for u1,x, and u2,x are also included.

Table 2
ABCDEFGHIJ..
1tt0t0t0t0t1t1t1t1..tf
2xu1u2u1,xu2,xu1u2u1,xu2,x..u2,x
3xs
4x1
5x2u1(x2,t1)
....
Nxe

Transient View Format

In this format, the roles of x and t are interchanged. Output time points are displayed in the first column, whereas spatial points are displayed in the first row. This format allows you to easily plot transient views of the variables at desired spatial points.

This example is discussed in detail in the following publication:
Ghaddar C.K., Rapid Modeling and Parameter Estimation of Partial Differential Algebraic Equations by a Functional Spreadsheet Paradigm.Math. Comput. Appl. 2018, 23, 39.

Consider the following coupled PDE equations which model a simplified battery system:

m11 Φ1t + m12 Φ2t = σx Φ1,xx - a x Φ1 - Φ2 - f t,x m21 Φ1t + m22 Φ2t + m23 Φ2t = κ x   Φ2,xx + a x Φ1 - Φ2 - f t,x m32 Φ2t + m33 Φ3t = σ x   Φ3,xx - a x ( Φ3 - Φ2 - f t,x

Where

Properties

σx=10xlc0lc<x<ls1lsxla,      κx=β*10xlc0lc<x<ls1lsxla,      ax=γ1060xlc0lc<x<ls106lsxla,      ft,x=3e-t0xlc0else

Parameters

lc=0.0002,      β=1,      ls=0.0002254,      γ=0.5,      la=0.0004254,      mc=106,      mij=106-mc0-mc106-mc0-mc106

Initial Values

Φ1(0,x)= Φ2(0,x)= Φ3(0,x)=0

Boundary conditions

x = 0 x = lc x = ls x = la
σ0 Φ1,x=1 Φ1,x=0 Φ3,x=0 Φ3=0
Φ2,x=0 κlcΦ2,x- = κlc Φ2,x+ κls Φ2,x-= κls Φ2,x+ Φ2,x=0

Solution

Working with named system variables, parameters, and property functions, we define the complete Excel model including equations, regions, and boundary conditions formulas as shown in Table 1.

Table 1
ABCDEFGH
1System Variables with initial valuesSystem Parameters
2tlc2.00E-04
3xls2.254E-04
4phi_10la4.254E-04
5phi_20beta1
6phi_30gamma0.5
7phi1xmc1E+06
8phi2x
9phi3xMass Matrix
10phi1xx1E+06=-mc0
11phi2xx=-mc1E+06=-mc
12phi3xx0=-mc1E+06
13
14Property Functionsx Domain0=la
15sigma=IF(x<=lc,1,IF(x<ls,0,1))t Interval01Boundary conditions matrix
16kappa=beta*IF(x<=lc,1,IF(x<ls,2,1))0R=sigma*phi1x-1
17a=gamma*IF(x<=lc,1E6,IF(x<ls,0,1E6))0N=phi2x
18f=IF(x<=lc,3*EXP(-2*t),0)=lcN=phi1x
19=lcC=kappa*phi2x
20System RHS EquationsRegions=lsC=kappa*phi2x
21eq1=sigma*phi1xx-a*(phi_1-phi_2-f)0=lc=lsN=phi3x
22eq2=kappa*phi2xx+a*(phi_1-phi_2-f)0=la=laN=phi2x
23eq3=sigma*phi3xx-a*(phi_3-phi_2-f)=ls=la=laD=phi_3

To aid readability, we assign assign names for selected Excel ranges in Table 1 which represent input arguments to PDASOLVE as shown in Table 2.

Table 2
Excel RangeAssigned Name
B21:B23Eqns
B2:B12Vars
F16:H23BCs
D14:E14xDom
D15:E15tDom
C10:E12M
C21:D23Rgns

Next, we evaluate the formula
=PDASOLVE(Eqns, Vars, BCs, xDom, tDom, M, Rgns, , {"FORMAT","TCOL1"})
in array G15:S27 and obtain the solution as shown in Table 3.

Table 3
GHIJKLMNOPQRS
15

x

0000.00020.00020.000200.00022540.00022540.00042540.00042540.00043
16

t

phi_1

phi_2

phi_3

phi_1

phi_2

phi_3

phi_1

phi_2

phi_3

phi_1

phi_2

phi_3

170000000000000
180.10.0490406-0.003126900.0491381-0.002664600-0.0026093-0.00039370-0.0022350
190.20.0888748-0.00478200.0889728-0.004394500-0.0043481-0.00033110-0.0040330
200.30.1204078-0.006139300.1205061-0.005812900-0.0057737-0.00028010-0.0055070
210.40.1451792-0.007253900.1452778-0.006977400-0.0069442-0.00023850-0.0067160
220.50.1644121-0.008169700.164511-0.007933900-0.0079055-0.00020460-0.007710
230.60.1791579-0.008924400.179257-0.008721700-0.0086973-0.00017690-0.0085280
240.70.1902228-0.009546800.190322-0.009371100-0.0093499-0.00015440-0.0092010
250.80.1983065-0.010062100.1984059-0.009908300-0.0098896-0.00013610-0.0097580
260.90.2039649-0.0104900.2040644-0.01035400-0.0103374-0.00012130-0.010220
2710.2076535-0.010846800.2077531-0.010725300-0.0107105-0.00010930-0.0106050

In Figure 1, we plot phi2(x,t) at the interface x = lc together with imperical values for Φ(lc,t) Clearly the default model does not agree well with experimental data. In Example 2, we will show you how to optimize the model parameters so it agrees with experimental data accurately.

Figure 1

Verification of Continuity Boundary Conditions

The result shown in Table 1, is insufficient to verify the flux continuity conditions κlcΦ2,x- = κlc Φ2,x+ and κls Φ2,x-= κls Φ2,x+ are staisfied at the regions interfaces x = lc and x = ls. Here, we illustrate use of additional optional parameters to report the 1st derivatives at custom spatial points just before and after x = lc and x = ls, as well as, enhance the solution accuracy.

Table 4 shows the definition of new optional parameters consisting of: custom output points in G20:N20 which is named as xOut; relative tolerances in D21:D23 which is named as rTols; and a set of key/value pairs in A21:B23 which is named Cntrls. The keys include NDRVOUT with a value of 1, to report the first derivative in the output, and MAX_GRID_SPACING which ensures that the maximum distance between computational spatial grid nodes does not exceed the specified value of 1.0e-5.

Table 4
ABCDEFGHIJKLMN
20Solver OptionsRel. Tolerancexout0=lc-0.001*lc=lc=lc+0.001*lc=ls-0.001*ls=ls=ls+0.001*ls=la
21FORMATTCOL10.0001
22NDRVOUT10.000001
23MAX_GRID_SPACING0.000010.000001

We evaluate the enhanced formula
=PDASOLVE(Eqns, Vars, BCs, xOut, tDom, M, Rgns, rTols, Cntrls)
in array A50:AW62 to obtain a new solution. We show in Table 5 values for phi2x just before and after the interfaces at x = lc and x = ls. A quick inspection of the ratios Φ2,x-/Φ2,x+ = κlc+/κlc-=2, and Φ2,x-/Φ2,x+ = κls+/κls-=0.5 in Table 5 shows that the flux continuity conditions are satisfied at all times

Table 5
LRXADAJAP
50x0.00019980.00020.0002002..0.000225170.0002250.000225625
51tphi2Xphi2xphi2xphi2xphi2xphi2x
5200Undef00Undef0
530.10.999180576Undef0.499576930.44674817Undef0.891514551
540.20.999180576Undef0.499595290.44905949Undef0.89617388
550.30.999180577Undef0.499612860.45127134Undef0.900632825
560.40.999180568Undef0.499629670.45338715Undef0.904898145
570.50.999180579Undef0.499645750.45541108Undef0.908978265
580.60.999180562Undef0.499661130.45734714Undef0.912881228
590.70.999180583Undef0.499675850.45919913Undef0.916614717
600.80.999180566Undef0.499689920.46097071Undef0.920186109
610.90.999180582Undef0.499703390.46266537Undef0.92360242
6210.999180573Undef0.499716270.46428644Undef0.926870399

This example is a continuation of Example 1.

The computed values for phi2 at the interface x = lc based on the default system parameters, are clearly in disagreement with the observed data as shown in Figure 1 of Example 1. The observed values for Φ(lc,t) are give in Table 6.

Table 6
MN
1timeEmperical phi2
200
30.1-0.002052705
40.2-0.00380538
50.3-0.005944579
60.4-0.006668727
70.5-0.008436437
80.6-0.009136398
90.7-0.009734491
100.8-0.010031262
110.9-0.010050529
121-0.010387955

Our task is to compute optimal values for the design parameters γ, and mc which minimize the sum of square residuals between our model predictions and the observed data in Table 6. We have two options to solve this dynamical minimization problem: Excel built-in NLP Solver; and NLSOLVE. We demonstrate both procedures below.

Using Excel Solver

We define an objective formula to be minimized. The objective formula calculates the summation of the square residuals between the observed values and their corresponding computed values from the solution result in Table 3 of Example 1. The definition of the objective formula is shown in S2 of Table 7 with initial value of 0.00033.

Table 7
PS
1=(L18-N2)^2Objective
2=(L19-N3)^2=SUM(P1:P10)
3=(L20-N4)^2
4=(L21-N5)^2
5=(L22-N6)^2
6=(L23-N7)^2
7=(L24-N8)^2
8=(L25-N9)^2
9=(L26-N10)^2
10=(L27-N11)^2
10=(L28-N12)^2

The remaining task is to configure and run Excel Solver command. We configure the solver exactly as shown in Figure 2. We select to minimize the objective formula S2, by varying gamma and mc, subject to the upper bound constraint mc <= 10e6 which is added for numerical stability since larger values for the off-diagonal mass matrix element mc leads to a divergent system.

Figure 2

We run the solver which reports a feasible solution in less than a minute of computing time as shown in the Answer Report generated by the Excel Solver in Figure 3. Upon accepting the solution, Excel automatically updates the spreadsheet to reflect the optimal computed values. The optimized solution for phi2 is shown in Figure 4 exhibiting excellent agreement with observed values.

Figure 3
Figure 4

Using NLSOLVE

Unlike Excel Solver which requires an objective formula, NLSOLVE requires the list of individual residual constraint formulas which we define as shown in Table 8.

Table 8
R
1=ARRAYVAL(L18)-N2
2=ARRAYVAL(L19)-N3
3=ARRAYVAL(L20)-N4
4=ARRAYVAL(L21)-N5
5=ARRAYVAL(L22)-N6
6=ARRAYVAL(L23)-N7
7=ARRAYVAL(L24)-N8
8=ARRAYVAL(L25)-N9
9=ARRAYVAL(L26)-N10
10=ARRAYVAL(L27)-N11
10=ARRAYVAL(L27)-N12

What is ARRAYVAL? ARRAYVAL is one of the criterion functions that allow us to define dynamic constraints on the solution array for input to NLSOLVE. In this case we are simply picking a value from the solution but in general, a criterion function can be used to specify more complex constraint on the solution.

Next, we evaluate the formula =NLSOLVE(R1:R10,(gamma,mc)) in array S4:T8 and obtain virtually the same solution found by Excel Solver as shown in Table 9. However, NLSOLVE, converges in fewer iterations in less than one third of the time. Since NLSOLVE is a pure spreadsheet function, it does not modify its arguments or update the spreadsheet calculation. Therefore, to update the spreadsheet result, the computed values for gamma and mc must copied manually to their respective variable cells D6 and D7 of Table 3 of Example 1.

Table 9
ST
4gamma0.195638747
5mc958050.8287
6SSERROR1.61147E-06
7ITRN10
8TIME (s)18.729

PDASOLVE implements the method of lines. Spatial discretization is carried out using a second-order variable-step finite difference scheme, and the resulting implicit DAE system is integrated by RADAU5 with adaptive step control.

  • Ghaddar C.K., Rapid Modeling and Parameter Estimation of Partial Differential Algebraic Equations by a Functional Spreadsheet Paradigm.Math. Comput. Appl. 2018, 23, 39.
  • Schiesser W.E. The Numerical Method of Lines, San Diego, CA: Academic Press, 1991.
  • Ernst Hairer and Gerhard Wanner, Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, Springer Series in Computational Mathematics, 1996.
  • J.W. Baish, K. Mukundakrishnan and P.S. Ayyaswamy, Numerical Models of Blood Flow Effects in Biological Tissues, Advances in Numerical Heat Transfer Volume 3: Numerical Implementation of Bioheat Transfer Models and Equations, CRC Press, Taylor and Francis Group, W.J. Minkowycz and E.M. Sparrow, Eds., Chapter 2, p. 29-74, 2009.
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