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:

Where

Properties

$\sigma \left(x\right)=\left\{\begin{array}{cc}1& 0\le x\le l_c\\ 0& l_c<x<l_s\\ 1& l_s\le x\le l_a\end{array}\right.$,      $\kappa \left(x\right)=\beta *\left\{\begin{array}{cc}1& 0\le x\le l_c\\ 0& l_c<x<l_s\\ 1& l_s\le x\le l_a\end{array}\right.$,      $a\left(x\right)=\gamma \left\{\begin{array}{cc}{10}^6& 0\le x\le l_c\\ 0& l_c<x<l_s\\ {10}^6& l_s\le x\le l_a\end{array}\right.$,      $f\left(t,x\right)=\left\{\begin{array}{cc}3e^{-t}& 0\le x\le l_c\\ 0& else\end{array}\right.$

Parameters

Initial Values

${\Phi }_1(0,x)=$ ${\Phi }_2(0,x)=$ ${\Phi }_3(0,x)=0$

Boundary conditions

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.

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.

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.

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 = l_c and x = l_s. Here, we illustrate use of additional optional parameters to report the 1st derivatives at custom spatial points just before and after x = l_c and x = l_s, 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.

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 = l_c and x = l_s. 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

This example is a continuation of Example 1.

The computed values for phi2 at the interface x = l_c 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 Φ(l_c,t) are give in Table 6.

Our task is to compute optimal values for the design parameters γ, and m_c 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.

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 m_c leads to a divergent system.

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.

---

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.

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.