REPRESENTATION OF SOLUTIONS TO LINEAR ALGEBRAIC EQUATIONS

by Fred Krylov

A dissertation submitted to the faculty of

The University of Utah

in partial fulfillment of the requirements for the degree of

Doctor of Philosophy

Department of Mathematics

The University of Utah March 1993 Copyright \begin{matrix}c \\ ◯ \end{matrix} Fred Krylov 1993

All Rights Reserved

THE UNIVERSITY OF UTAH GRADUATE SCHOOL

SUPERVISORY COMMITTEE APPROVAL

of a dissertation submitted by

Fred Krylov

This dissertation has been read by each member of the following supervisory committee and by majority vote has been found to be satisfactory.

Chair: Fletcher Gross Hans Othmer Jim Carlson Grant Gustafson Nick Korevaar THE UNIVERSITY OF UTAH GRADUATE SCHOOL

FINAL READING APPROVAL

To the Graduate Council of the University of Utah:

I have read the dissertation of  \underline{Fred Krylov}  in its final form

and have found that (1) its format, citations, and bibliographic style are consistent and acceptable; (2) its illustrative materials including figures, tables, and charts are in place; and (3) the final manuscript is satisfactory to the Supervisory Committee and is ready for submission to The Graduate School.

Date Fletcher Gross

Chair, Supervisory Committee

Approved for the Major Department

Paul Fife

Chair/Dean

Approved for the Graduate Council

Ann W. Hart

Dean of The Graduate School

ABSTRACT

This dissertation is described here in less than 350 words using no footnotes, diagrams, references or outside anything.

For my cat, Mouse, a few lines only.

CONTENTS LIST OF FIGURES LIST OF TABLES

NOTATION AND SYMBOLS

Most of the following may be found in [?].

n  n-dimensional Euclidean space.
Ω  A bounded open subset of ^{n}.
∂Ω  The boundary of Ω.
B_{r}  = {x : |x| < r}, the ball of radius r.
D^{β}f  = ∂^{|β|}f/∂x_{β_{1}}∂x_{β_{2}} ⋯ ∂x_{β_{n}}, |β| ≡ ∑_{i} β_{i}.
∇f  = (∂f/∂x_{1}, ∂f/∂x_{2}, …, ∂f/∂x_{n}), the gradient of f.
div{g}  = ∂g_{1}/∂x_{1} + ∂g_{2}/∂x_{2} + … + ∂g_{n}/∂x_{n}, the divergence of g.
Δf  = div{∇f}, the Laplacian of f.
C^{k}(Ω)  Functions defined on Ω which have k continuous derivatives.
C_{0}^{k}(Ω)  C^{k}(Ω) functions which vanish at the boundary.
C^{0,α}(Ω)  Hölder continuous functions with Hölder constant α.
C^{k,α}(Ω)  C^{k}(Ω) functions with C^{0,α}(Ω) derivatives (up to order k).
C_{0}^{k,α}(Ω)  C^{k,α}(Ω) functions which vanish at the boundary.
∥f∥_{L^{p}(Ω)}  = (∫_{Ω} |f|^{p}dx)^{1/p}, the L^{p} norm.
L^{p}(Ω)  The space of p integrable functions (the L^{p} norm is bounded).
∥f∥_{W^{k,p}(Ω)} = (∑_{|β|≤k} ∫_{Ω} |D^{β}f|^{p}dx)^{1/p}, the Sobolev norm.
W^{k,p}(Ω)  The space of functions with bounded Sobolev norm.
W_{0}^{k,p}(Ω)  W^{k,p}(Ω) functions that vanish a.e. at the boundary.

ACKNOWLEDGEMENTS

This page is optional. It’s in the table of contents and it’s labeled ’ACKNOWLEDGMENTS’ even though the spelling ‘ACKNOWLEDGEMENTS’ is also correct. This page should be at the end of the preface, if one exists, or a separate page, if no preface is used.

Recognize faculty, friends, family but not the typist, editor or persons paid to produce the manuscript from draft copy. Recognize funding, special permission to reproduce figures and copyrighted material.

CHAPTER 1

INTRODUCTION

1.1 The Sample.tex file

1.1.1 Ordinary Text

The ends of words and sentences are marked by spaces. It doesn’t matter how many spaces you type; one is as good as 100. The end of a line counts as a space.^{1}

One or more blank lines denote the end of a paragraph.

Since any number of consecutive spaces are treated like a single one, the formatting of the input file makes no difference to LaTEX, but it makes a difference to you. When you use LaTEX, making your input file as easy to read as possible will be a great help as you write your document and when you change it. This sample file shows how you can add comments to your own input file.

Because printing is different from typewriting, there are a number of things that you have to do differently when preparing an input file than if you were just typing the document directly. Quotation marks like “this” have to be handled specially, as do quotes within quotes: “ ‘this’ is what I just wrote, not ‘that’ .”

Dashes come in three sizes: an intra-word dash, a medium dash for number ranges like 1–2, and a punctuation dash—like this.

A sentence-ending space should be larger than the space between words within a sentence. You sometimes have to type special commands in conjunction with punctuation characters to get this right, as in the following sentence. Gnats, gnus, etc. all begin with G. You should check the spaces after periods when reading your output to make sure you haven’t forgotten any special cases. Generating an ellipsis … with the right spacing around the periods requires a special command.

LaTEX interprets some common characters as commands, so you must type special commands to generate them. These characters include the following: $ & % # { and }.

1This is a sample input file. Comparing it with the output it generates can show you how to produce a simple document of your own.

In printing, text is emphasized by using an italic type style. A long segment of text can also be emphasized in this way. Text within such a segment

given additional emphasis with Roman type. Italic type loses its ability to emphasize and become simply distracting when used excessively.

It is sometimes necessary to prevent LaTEX from breaking a line where it might otherwise do so. This may be at a space, as between the “Mr.” and “Jones” in “Mr. Jones,” or within a word—especially when the word is a symbol like itemnum that makes little sense when hyphenated across lines.

Footnotes^{2} pose no problem. LaTEX is good at typesetting mathematical formulas like x − 3y = 7 or

a_{1} > x^{2n}/y^{2n} > x′.

Remember that a letter like x is a formula when it denotes a mathematical symbol, and should be treated as one.

1.1.2 Displayed Text

Text is displayed by indenting it from the left margin.

1.1.2.1 Quotations

Quotations are commonly displayed. There are short quotations

This is a short a quotation. It consists of a single paragraph of text. There is no paragraph indentation.

and longer ones.

This is a longer quotation. It consists of two paragraphs of text. The beginning of each paragraph is indicated by an extra indentation.

This is the second paragraph of the quotation. It is just as dull as the first paragraph.

1.1.2.2 Lists

Another frequently-displayed structure is a list.

2This is an example of a footnote.

1.1.2.2.1 Itemize. The following is an example of an itemized list.

• This is the first item of an itemized list. Each item in the list is marked with a “tick”. The document style determines what kind of tick mark is used.

• This is the second item of the list. It contains another list nested inside it. The inner list is an enumerated list.

1. This is the first item of an enumerated list that is nested within the itemized list.

2. This is the second item of the inner list. LaTEX allows you to nest lists deeper than you really should.

This is the rest of the second item of the outer list. It is no more interesting than any other part of the item.

• This is the third item of the list.

1.1.2.2.2 Verse. You can even display poetry.

There is an environment for verse Whose features some poets will curse. For instead of making Them do all line breaking, It allows them to put too many words on a line when they’d

rather be forced to be terse.

1.1.2.3 Mathematics

Mathematical formulas may also be displayed. A displayed formula is one-line long; multiline formulas require special formatting instructions.

x′ + y^{2} = z_{i}^{2}

Don’t start a paragraph with a displayed equation, nor make one a paragraph by itself.

1.2 More examples: Jeff McGough’s Thesis

Equations like γ = 0 that don’t need numbering may be set inline by the coding $\gamma = 0$ or displayed by

$$
\gamma = 0.
$$

Numbered equations are set as shown in the next paragraph. They use the theorem

environments defined in thesis.sty:

\newtheorem{thrm}{Theorem}
\newtheorem{lem}[thrm]{Lemma}
\newtheorem{cor}[thrm]{Corollary}
\newtheorem{rem}[thrm]{Remark}
\newtheorem{defn}[thrm]{Definition}
\newtheorem{exmpl}[thrm]{Example}

The Gelfand problem is the following elliptic boundary value problem:

\begin{matrix}Δu + λe^{u} = 0, & u ∈ Ω, \\ u = 0, & u ∈ ∂Ω. \end{matrix}  (1.1)

The previous equation had a label. It may be referenced as equation (??).

1.3 History of the Gelfand problem

According to Bebernes and Eberly [?, p.46], Gelfand was “the first to make an in-depth

study” of (??). Following this statement they briefly outline the history of the Gelfand

problem.

For dimension n = 1, Liouville [?] first studied and found an explicit solution in 1853. For n = 2, Bratu [?] found an explicit solution in 1914. Frank-Kamenetski [?] rediscovered these results in his development of thermal explosion theory. Joseph and Lundgren [?] gave an elementary proof via phase plane analysis of the multiple existence of solutions for dimensions n ≥ 3.

From Zeidler [?]:

\begin{aligned}div j & = & f, & x ∈ Ω, \\ u & = & g_{1}, & x ∈ ∂Ω_{1}, \\ jν & = & g_{2}, & x ∈ ∂Ω_{2}, \end{aligned}  (1.2)

where

j = h(|∇u|^{2})∇u  (1.3)

and Ω is a bounded domain in ^{n} with smooth boundary ∂Ω = \overline{∂Ω_{1}} ∪\overline{∂Ω_{2}}, ∂Ω_{1} ∩∂Ω_{2} = ∅

and ν is the normal vector to ∂Ω.

Lemma 1 Assuming that ∂Ω_{2} = ∅ and that h(t) = 1, we have

\begin{matrix}Δu = f, & x ∈ Ω, \\ u = g_{1}, & x ∈ ∂Ω. \end{matrix}

Corollary 2 If g_{2} = 0 then

\begin{matrix}Δu = f, & x ∈ Ω, \\ u = 0, & x ∈ ∂Ω. \end{matrix}

1.4 Fundamental results

The investigation of the Gelfand problem begins with examining the ..... (this paragraph continues for many lines).

Theorem 3 (Joseph-Lundgren [?]) Boundary value problem (??) has positive radial solutions u on the unit ball which depend on n and λ in the following manner.

1. For n = 1, 2, there exists λ^{∗} > 0 such that (a) for 0 < λ < λ^{∗} there are two positive solutions,

(b) for λ = λ^{∗} there is a unique solution, and

(c) for λ > λ^{∗} there are no solutions.

2. For 3 ≤ n ≤ 9, let \overline{λ} = 2(n − 2); then there exist positive constants λ_{∗}, λ^{∗} with 0 < λ_{∗} < \overline{λ} < λ^{∗}, such that

(a) for λ = λ^{∗} there is a unique solution,

(b) for λ > λ^{∗} there are no solutions,

(c) for λ = \overline{λ} there is a countably infinite number of solutions,

(d) for λ ∈ (λ_{∗}, λ^{∗}), λ ≠ \overline{λ}, there is a finite number of solutions,

(e) for λ < λ_{∗} there is a unique solution.

3. For n ≥ 10, let λ^{∗} = 2(n − 2) then (a) for λ ≥ λ^{∗} there are no solutions,

(b) for λ ∈ (0, λ^{∗}) there is a unique solution.

QUADRATIC NONLINEARITIES

In this chapter we derive results for the quadratic equation.

2.1 Derivation of the quadratic formula

A quadratic equation is one of the form

ax^{2} + bx + c = 0  (2.1)

where a, b, c are known constants and x is the unknown. The results are summarized in Table ?? and Table ?? below.

2.2 Application of the quadratic formula

If the differential operator generates a nonnegative form, then an inequality is based on the following considerations. See Figure ?? for n = 1, 2, Figure ?? for 3 ≤ n ≤ 9 and Figure ?? for n ≥ 10.

Table 2.1. PDE solve times, 15^{3} + 1 equations.

Precond. Time  Nonlinear Krylov  Function Precond.
Iterations Iterations calls  solves
None 1260.9u 3 26 30 0

(21:09)

FFT 983.4u 2 5 8 7

(16:31)

MILU 629.7u 3 11 15 14

(10:36)

Table 2.2. Convergence properties of RQI.

Object  Normal Matrices Diagonalizable Matrices Defective Matrices
ρ  Stationary at ev’s. Stationary at ev’s.  Stationary at ev’s.
∥r_{k}∥  → 0 as k → ∞.  Can oscillate.  Can oscillate.
ρ_{k}  Converges.  Unknown.  Unknown.
Convergence to is cubic.  is quadratic.  is linear.

eigensets

∥u∥  λ

Figure 2.1. Gelfand equation on the ball, n = 1, 2.

∥u∥  λ

Figure 2.2. Gelfand equation on the ball, 3 ≤ n ≤ 9.

∥u∥  λ

Figure 2.3. Gelfand equation on the ball, n ≥ 10.

SYSTEMS

3.1 Diagrams made with diagram.sty

An example diagram appears below in Figure ??. This is typical of what can made with the diagram package.

3.2 Sample diagrams from diagram.tex

Example diagrams reproduced here were taken from various sources. Compare the three diagrams of increasing sizes in Figure ??, Figure ??, Figure ?? with the three diagrams in Figure ??, Figure ??, Figure ??.

[figure: U; X; Y − ∂Q; Y; i_{1}; π; j_{1}]

Figure 3.1. Diagram example

[figure: A; B; C; D]

Figure 3.2. Base diagram, Arrowlength = 3.0em

[figure: A; B; C; D]

Figure 3.3. Same as Figure ??, but Arrowlength = 6.0em

[figure: A; B; C; D]

Figure 3.4. Same as Figure ??, but Arrowlength = 12.0em

[figure: A; B; C; D]

Figure 3.5. Base figure, same as Figure ??.

[figure: A; B; C; D]

Figure 3.6. Same as Figure ??, but Bignode = strut hspace 6.0em.

[figure: A; B; C; D]

Figure 3.7. Same as Figure ??, but Bignode = strut hspace 12.0em

Below we show diagrams from the manual with a few modifications. The first in Figure ?? is essentially as it appears in the manual, whereas the second, Figure ?? has been rescaled to a larger size.

Below are several diagrams created by Bill Richter. The first, Figure ?? is modified slightly to produce Figure ??. Both use fractur fonts. The last one, Figure ??, is a complicated example illustrating the limits of what can be done with diagrams.

The diagram below in Figure ??, the last of our series of illustrations, is by Anders Thorup (thorup@math.ku.dk), originally done with a package developed by himself and

Steven Kleiman (kleiman@math.mit.edu):

[figure: A; B^{∗}; C; D; E; a; c; u; b^{∗}; d; v; e]

Figure 3.8. First diagram from manual

[figure: A; B^{∗}; C; D; E; a; c; u; b^{∗}; d; v; e]

Figure 3.9. First diagram from manual, rescaled.

[figure: K ( /2, 8n + 1); K ( /4, 8n − 1); E; ΣΩX ∧ ΩX; ΣΩX; X; K ( /4, 8n).; Sq^{2}; Θ; π; H_{μ}; σ(α⊗α); σ; σ̃; α^{2}]

Figure 3.10. Bill Richter, first diagram

[figure: ΩΣA; Ω^{2} (ΣA ∧ ΣA); #; ΩB; F; Ω^{2} (B ∧ ΣA); A; X; B.; λ_{2}; δ; Ω(∂); 𝔥; π; J; Ω^{2}(∂∧id); 𝔢; f; E; h]

Figure 3.11. Bill Richter, second diagram

[figure: ΩS^{5}; ; Ω (M_{2⊂}^{5}); ΩS^{5}; G; J (S^{4} ∧ S^{4}); J_{2} (M_{2⊂}^{4}); G_{2}; S^{8}; ℵ; M_{2⊂}^{12}; S^{4}; M_{2⊂}^{5}; S^{5}; Ω(pinch); δ; Ω(2⊂); 𝔥_{2}; π; J; δ_{2}; ⊂; 𝔥_{2}; ⊂; E; τ̄; τ; ⊂; 𝔢; E; pinch]

Figure 3.12. Bill Richter, third diagram

[figure: H^{k}(B_{G} × N; Q) = H_{G}^{k}(N; Q); H_{G}^{k}(F_{j}; Q); H_{G}^{k}(M; Q); H^{k}(N; Q); H^{k}(F_{j}; Q); H^{k}(M; Q); f_{j}^{∗}; p^{∗}; f̃^{∗}; q_{j}^{∗}; i_{j}^{∗}; i^{∗}; f̃_{j}^{∗}=f_{j}^{∗}; f̃^{∗}=f^{∗}; i_{j}^{∗}]

Figure 3.13. Anders Thorup diagram

APPENDIX

CLASSICAL IDENTITIES

Rellich’s identity

Standard developments of Pohozaev’s identity used an identity by Rellich [?], reproduced here.

Lemma 1 (Rellich) Given L in divergence form and a, d defined above, u ∈ C^{2}(Ω), we have

∫_{Ω}(−Lu)∇u · (x − \overline{x}) dx = (1 − \frac{n}{2}) ∫_{Ω} a(∇u, ∇u) dx − \frac{1}{2} ∫_{Ω} d(∇u, ∇u) dx  (A.1)

+\frac{1}{2} ∫_{∂Ω} a(∇u, ∇u)(x − \overline{x}) · ν dS − ∫_{∂Ω} a(∇u, ν)∇u · (x − \overline{x}) dS.

Proof: There is no loss in generality to take \overline{x} = 0. First rewrite L:

Lu = \frac{1}{2} [∑_{i} ∑_{j} \frac{∂}{∂x_{i}} (a_{ij} \frac{∂u}{∂x_{j}} ) + ∑_{i} ∑_{j} \frac{∂}{∂x_{i}} (a_{ij} \frac{∂u}{∂x_{j}} )]

Switching the order of summation on the second term and relabeling subscripts, j → i and i → j, then using the fact that a_{ij}(x) is a symmetric matrix, gives the symmetric form needed to derive Rellich’s identity.

Lu = \frac{1}{2} ∑_{i,j} [\frac{∂}{∂x_{i}} (a_{ij} \frac{∂u}{∂x_{j}} ) + \frac{∂}{∂x_{j}} (a_{ij} \frac{∂u}{∂x_{i}} )] .  (A.2)

Multiplying −Lu by \frac{∂u}{∂x_{k}} x_{k} and integrating over Ω, yields

∫_{Ω}(−Lu)\frac{∂u}{∂x_{k}} x_{k} dx = −\frac{1}{2} ∫_{Ω} ∑_{i,j} [\frac{∂}{∂x_{i}} (a_{ij} \frac{∂u}{∂x_{j}} ) + \frac{∂}{∂x_{j}} (a_{ij} \frac{∂u}{∂x_{i}} )] \frac{∂u}{∂x_{k}} x_{k} dx

Integrating by parts (for integral theorems see [?, p. 20]) gives

= \frac{1}{2} ∫_{Ω} ∑_{i,j} a_{ij} [\frac{∂u}{∂x_{j}} \frac{∂^{2}u}{∂x_{k}∂x_{i}} + \frac{∂u}{∂x_{i}} \frac{∂^{2}u}{∂x_{k}∂x_{j}} ] x_{k} dx

+\frac{1}{2} ∫_{Ω} ∑_{i,j} a_{ij} [\frac{∂u}{∂x_{j}} δ_{ik} + \frac{∂u}{∂x_{i}} δ_{jk}] \frac{∂u}{∂x_{k}} dx

−\frac{1}{2} ∫_{∂Ω} ∑_{i,j} a_{ij} [\frac{∂u}{∂x_{j}} ν_{i} + \frac{∂u}{∂x_{i}} ν_{j}] \frac{∂u}{∂x_{k}} x_{k} dx

= I_{1} + I_{2} + I_{3}, where the unit normal vector is ν. One may rewrite I_{1} as

I_{1} = \frac{1}{2} ∫_{Ω} ∑_{i,j} a_{ij} \frac{∂}{∂x_{k}} (\frac{∂u}{∂x_{i}} \frac{∂u}{∂x_{j}} ) x_{k} dx

Integrating the first term by parts again yields

I_{1} = −\frac{1}{2} ∫_{Ω} ∑_{i,j} a_{ij} (\frac{∂u}{∂x_{i}} \frac{∂u}{∂x_{j}} ) dx + \frac{1}{2} ∫_{∂Ω} ∑_{i,j} a_{ij} (\frac{∂u}{∂x_{i}} \frac{∂u}{∂x_{j}} ) x_{k}ν_{k} dS

−\frac{1}{2} ∫_{Ω} ∑_{i,j} (\frac{∂u}{∂x_{i}} \frac{∂u}{∂x_{j}} ) x_{k} \frac{∂a_{ij}}{∂x_{k}} dx.

Summing over k gives

∫_{Ω}(−Lu)(∇u · x) dx = −\frac{n}{2} ∫_{Ω} ∑_{i,j} a_{ij} (\frac{∂u}{∂x_{i}} \frac{∂u}{∂x_{j}} ) dx

+\frac{1}{2} ∫_{∂Ω} ∑_{i,j} a_{ij} (\frac{∂u}{∂x_{i}} \frac{∂u}{∂x_{j}} ) (x · ν) dS − \frac{1}{2} ∫_{Ω} ∑_{i,j} (\frac{∂u}{∂x_{i}} \frac{∂u}{∂x_{j}} ) (x · ∇a_{ij}) dx

+\frac{1}{2} ∫_{Ω} ∑_{i,j,k} a_{ij} [\frac{∂u}{∂x_{j}} \frac{∂u}{∂x_{k}} δ_{ik} + \frac{∂u}{∂x_{i}} \frac{∂u}{∂x_{k}} δ_{jk}] dx

−\frac{1}{2} ∫_{∂Ω} ∑_{i,j} a_{ij} [\frac{∂u}{∂x_{j}} ν_{i} + \frac{∂u}{∂x_{i}} ν_{j}] (∇u · x) dS.

Combining the first and fourth term on the right-hand side simplifies the expression

∫_{Ω}(−Lu)(∇u · x) dx = (1 − \frac{n}{2}) ∫_{Ω} ∑_{i,j} a_{ij} (\frac{∂u}{∂x_{i}} \frac{∂u}{∂x_{j}} ) dx

+\frac{1}{2} ∫_{∂Ω} ∑_{i,j} a_{ij} (\frac{∂u}{∂x_{i}} \frac{∂u}{∂x_{j}} ) (x · ν) dS − \frac{1}{2} ∫_{Ω} ∑_{i,j} (\frac{∂u}{∂x_{i}} \frac{∂u}{∂x_{j}} ) (x · ∇a_{ij}) dx

−\frac{1}{2} ∫_{∂Ω} ∑_{i,j} a_{ij} [\frac{∂u}{∂x_{j}} ν_{i} + \frac{∂u}{∂x_{i}} ν_{j}] (∇u · x) dS.

Using the notation defined above, the result follows.

Fortran code

subroutine rhs(neq,v,rhsf) save

c c This subroutine computes the function values. Inputs are neq and c v, and on output the values of f are stored in the array of rhsf c

include 'parabolic.inc'

integer neq  10

integer i integer j integer k integer ind integer inde integer indw integer indn integer inds integer ind0

integer ind1  20

integer ind2

double precision v(neq) double precision rhsf(neq) double precision u(nv) double precision diff double precision diffn double precision diffxn double precision diffyn

double precision nl  30

c  write(*,*)'funct begin'

c c Compute F for the local dynamics, written as F(u)= −du/dt + f(u) c c c the system parameters c

c  p1  ! parameter F  40
c  p2  ! parameter k

do j = 1, ny

do i = 1, nx

c c set up index c

ind = (i−1)*nv + (j−1)*meq

c

c Extract the jth component at current time  50

nl = v(1+ind)*v(2+ind)*v(2+ind)

rhsf(1+ind) = (− nl + p1*(1.0d0 − v(1+ind)))*local

rhsf(2+ind) = ( nl − (p1+p2)*v(2+ind))*local

end do end do

c −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−  60

c c add diffusion for all species (zero diffusion c coefficient takes care of those that do not diffuse). c c −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−

do j = 1, ny do i = 1, nx

c indexing  70

c ind0 = (i−1)*nv + (j−1)*meq ! point indw = (i−2)*nv + (j−1)*meq ! west point inde = (i)*nv + (j−1)*meq ! east point indn = (i−1)*nv + (j)*meq ! north point inds = (i−1)*nv + (j−2)*meq ! south point

if(i.eq.1) indw = (nx−1)*nv + (j−1)*meq

if(i.eq.nx) inde = (j−1)*meq

if(j.eq.1) inds = (i−1)*nv + (ny−1)*meq  80

if(j.eq.ny) indn = (i−1)*nv

do k = 1, 2

c c First compute the contribution within a row at the current time c and at the preceding time. c

ind = k + ind0

ind1 = k + indw  90

ind2 = k + inde

diffxn = v(ind1) − 2.0d0*v(ind) + v(ind2)

c c Compute the contribution from the columns c

ind1 = k + indn ind2 = k + inds

100

diffyn = v(ind1) − 2.0d0*v(ind) + v(ind2) c Multiply by other factors and sum c

diff = d(k)*hxx*(diffxn + diffyn)*diffus

rhsf(ind) = rhsf(ind) + diff

110

end do end do

end do

return end