FEATURES OF IMPLEMENTING THE METHOD OF LINES ON AN ANALOG COMPUTER

I. G. Vintizenko

(Presented by the scientific seminar of the computing laboratory)

In modeling partial differential equations on analog computers (ACs), the most convenient method is the method of lines (Eterman, 1957). However, to date this method has not been investigated for all types of equations and boundary conditions, and, in solving specific problems, it may not only lead to significant errors but also completely distort the solution (Giloi, 1963).

The present work shows a feature of implementing the method of lines for the equation of nonstationary diffusion in a one-dimensional region with concentrations of the diffusing substances \(C_R\) and \(C_0\).

For zone 1 the diffusion equation (cylindrical coordinate system) is

\[ \frac{\partial C_R(r,t)}{\partial t} = D_R\left[ \frac{\partial^2 C_R(r,t)}{\partial r^2} + \frac{2}{r_0}\frac{\partial C_R(r,t)}{\partial r} \right], \tag{1} \]

\[ 0 \leq r \leq r_0 . \]

For zone 2

\[ \frac{\partial C_0(r,t)}{\partial t} = D_0\left[ \frac{\partial^2 C_0(r,t)}{\partial r^2} + \frac{2}{r_0^*}\cdot\frac{\partial C_0(r,t)}{\partial r} \right], \tag{2} \]

\[ r_0^* \leq r \leq \infty, \]

where the boundaries of the zones \(r_0\) and \(r_0^*\) may either coincide \((r_0=r_0^*)\) or not coincide;

\(D_R\) and \(D_0\) are the corresponding diffusion coefficients.

Boundary conditions

\[ D_R\frac{\partial C_R(r,t)}{d r}\bigg|_{r=0}=0,\qquad D_0\frac{\partial C_0(r,t)}{\partial r}\bigg|_{r=\infty}=0. \tag{3} \]

Conditions for matching the solutions at the boundary of the zones

\[ D_R\frac{\partial C_R(r,t)}{\partial r}\bigg|_{r=r_0}=q(t),\qquad D_0\frac{\partial C_0(r,t)}{\partial r}\bigg|_{r=r_0^*}=-q(t). \tag{4} \]

Initial conditions

\[ C_R(r,0)=\operatorname{const}_1,\quad C_0(r,0)=\operatorname{const}_2. \tag{5} \]

Additional condition imposed on the matching of the solutions,

\[ \theta C_R(t,r_0)=\exp(-wt)C_0(t,r_0^*), \tag{6} \]

where \(\theta\) and \(w\) are constants.

We divide the solution region into finite intervals with respect to the variable \(r\) (Fig. 1). For zone 1 \((0 \leq r \leq r_0)\) we take a constant step

Fig. 1. Division of the solution region into finite intervals

\[ h_1=\frac{r_0}{K+\frac{1}{2}}. \]

Then equation (1) is written in the form of a system of ordinary differential equations in total derivatives. For the \(i\)-th nodal point,

\[ \frac{dC_R^i(t)}{dt}+\gamma C_R^i(t)=\alpha C_R^{i+1}(t)+\beta C_R^{i-1}(t), \tag{7} \]

where

\[ \alpha=\frac{D_R}{r_0^2}\left(K+\frac{1}{2}\right)\left(K+\frac{3}{2}\right), \]

\[ \beta=\frac{D_R}{r_0^2}\left(K+\frac{1}{2}\right)\left(K-\frac{1}{2}\right), \]

\[ \gamma=\frac{2D_R}{r_0^2}\left(K+\frac{1}{2}\right)^2. \]

For the point \(i=0\), taking condition (3) into account,

\[ C_R^{-1}(t)=C_R^0(t), \]

\[ \frac{dC_R^0(t)}{dt}+\alpha C_R^0(t)=\alpha C_R^1(t). \]

In the region \(r_0^* \leq r \leq \infty\), replace the variable \(r\) by the variable \(\rho\) according to the formula:

\[ \rho=\frac{1}{r}. \]

\[ \frac{\partial}{\partial r}=-\rho^2\frac{\partial}{\partial \rho}, \]

\[ -\frac{\partial^2}{\partial r^2} = \rho^3\left[\rho\frac{\partial^2}{\partial \rho^2} +2\frac{\partial}{\partial \rho}\right], \]

in this case the partition step \(h_2\) will be equal to

\[ h_2=\frac{\rho_N}{N+\frac{1}{2}} = \frac{1}{r_0\left(N+\frac{1}{2}\right)}. \]

Then, for the \(j\)-th nodal point,

\[ \frac{dC_0^j(t)}{dt} = P C_0^{j-1}(t)-2Q_j C_0^j(t)+R_j C_0^{j+1}(t), \]

where

\[ Q_j=\left(j+\frac{1}{2}\right)\rho_j^2D_0\left(N+\frac{1}{2}\right), \]

\[ P_j= \left[ \frac{j+\frac{1}{2}}{N+\frac{1}{2}}(Q_j+2)-2 \right] \rho_j^2D_0\left(N+\frac{1}{2}\right), \]

\[ R_j= \left[ \frac{j+\frac{1}{2}}{N+\frac{1}{2}}(Q_j-2)+2 \right] \rho_j^2D_0\left(N+\frac{1}{2}\right). \]

For the point \(N=0\), since \(C_0^0(t)=C_0^{-1}(t)\) (condition 3),

\[ \frac{dC_0^0(t)}{dt} = (P-2Q_0)C_0^0(t)+R_0C_0^1(t). \]

The matching of the solutions \(C_R(t)\) and \(C_0(t)\) occurs at the point \(K=N\) (the counting of nodal points in both cases goes from the ends of the zones toward the common middle). The quantities \(q(t)\), \(C_R^K(t)\), and \(C_0^N(t)\) are found from equations (4) and (6), written in finite-difference form:

\[ \frac{D_R}{h_1}\left[C_R^{K-1}(t)-C_R^K(t)\right]=q(t), \]

\[ \frac{D_0}{r_0}\left(N+\frac{1}{2}\right) \left[C_0^{N-1}(t)-C_0^N(t)\right]=q(t), \]

\[ \theta C_R^K(t)=C_0^N(t)\exp(-wt). \]

The simulation circuits for the \(i\)-th and \(j\)-th nodes, as well as the circuit for obtaining \(q(t)\), are given in Fig. 2. During simulation, instability of the fluxes \(C^j(t)\) was detected at all nodal points; this also occurred in the stationary case when the matching circuit was disconnected (Fig. 2, c). The fluxes very rapidly increased to levels lying outside the linear range of the operational amplifiers.

A review of the literature showed that a similar effect had been noted (Giloi, 1963) in solving a simpler diffusion equation, but with the same boundary conditions. In (Giloi, 1963) it was proposed, in order to stabilize the scheme, to change the boundary conditions so that the gradient at the boundaries became nonzero.

In our case this method gave an incorrect solution, since although it was possible to stabilize the nodal potentials at the boundaries of the domain, the potentials of the intermediate points continued to increase, which did not correspond to the physical meaning of the problem.

Fig. 2. Simulation circuits:

a) circuit for simulating the \(j\)-th node of zone I of the solution domain;

b) circuit for simulating the \(j\)-th node of zone II of the solution domain;

c) circuit for simulating the gradient \(q(t)\) at the point \(K-N\).

A method was proposed for refining the solution, consisting in the fact that, in the steady state, the potentials of the \(i\)-th nodes were stabilized by decreasing the values of \(\alpha\) relative to the calculated ones. It turned out that quite definite values of \(\alpha\) at all \(i\) points corresponded to steady fluxes at all nodes. The experimentally found \(\alpha\) values differ from the calculated ones by 8%.

One of the variants of the nonstationary solution with refined values of \(\alpha\) coincided with one previously calculated on the Minsk-1 electronic digital computer; in this digital calculation, the partial differential equations were transformed into integral equations without any connection with the finite-difference approximation used in the simulation. The simulation was carried out on the EMU-10 analog computer, and the results were recorded by the ETsR-1 digital recorder.

Thus, a certain caution is required when modeling systems with distributed parameters by the method of lines. Setting the gradients equal to zero at both ends of the domain leads to instability of the finite-difference approximated system, which explicitly or implicitly distorts the solution. One way of combating this phenomenon is the experimental determination of a number of coefficients required to bring the system into a stationary state. Further investigations are needed to eliminate the peculiarities of the method of lines in solving problems on analog computers.

REFERENCES

I. I. Eterman. Continuous-Action Mathematical Machines. Mashgiz, Moscow, 1957.

Giloi W. Über stabilitätsschwierigkeiten bei der Lösung gewisser Fälle der Diffusions und Wärmeleitungsgleichung auf dem Analogrecher. (“On difficulties in obtaining, on analog computers, certain equations of diffusion and heat conduction.”)

Ann. assoc. internat. calcul. analog. 1963, 5, No. 2, pp. 95–98.