RHS of Laplace's Equation is f(u(x,y))

[beth92]

Homework Statement

I need to (computationally) solve the following linear elliptic problem for the function u(x,y):

\Delta u(x,y) = u_{x,x} + u_{y,y} = k u(x,y)

on the domain

\Omega = [0,1]\times[0,1] with u(x,y) = 1 at all points on the boundary.

Homework Equations

[/B]
I know that I have to use discretization to solve the problem by considering u(x,y) as an n x n grid of points separated by distance h (and mapping them to a one dimensional array of points to make a vector U of n² components). Then the relevant equation for each point U_i is:

(1) u_{x,x} + u_{y,y} = \frac{U_{i-1} + U_{i+1} +U_{i+n} + U_{i-n} - 4 U_{i}}{h^{2}} = k U_{i}

This is known to be true for this problem (see the link http://en.wikipedia.org/wiki/Discrete_Laplace_operator#Finite_Differences for the explanation)

and it results in a linear system of equations as follows:

(2) A U = kU + B

Where A is a sparse matrix which has five diagonals (1,1,-4,1,1) at the offsets ( -(n+1),-1,0,1,(n+1) ). Then B is another vector (of dimension n²) containing the boundary condition values - any point on the edge of the grid will contain a +1 term in equation (1) instead of the value of the point beside it (since U = 1 for all boundary points).

The system (2) should be solved by a linear solver (eg. conjugate gradient method or something similar) computationally. The result can be plotted as a surface and I believe all u(x,y) values should be between 0 and 1. (i.e. the surface should 'dip' down at all points within the boundaries.)

The Attempt at a Solution

[/B]
I understand the general theory behind writing the problem as a linear system of equations and that it should be solve-able as such. However when writing a program to solve this I'm unsure how to deal with

(a) The fact that the RHS of the linear system now depends on the solution U

(b) The boundary condition vector B

I'm used to solving linear systems of the form Ax=b where the right hand side is a constant, known vector. So my question is, how do I make the system (2) look like this? I will input a LHS and RHS but should all values of U be on the left, and which side does the vector B go on? Am I right in thinking that the RHS of the system must be constantly updated?

Any comments are appreciated, and I'm happy to clarify any other info about the problem. Thank you

 

[Ray Vickson]

Homework Statement

I need to (computationally) solve the following linear elliptic problem for the function u(x,y):

\Delta u(x,y) = u_{x,x} + u_{y,y} = k u(x,y)

on the domain

\Omega = [0,1]\times[0,1] with u(x,y) = 1 at all points on the boundary.

Homework Equations

[/B]
I know that I have to use discretization to solve the problem by considering u(x,y) as an n x n grid of points separated by distance h (and mapping them to a one dimensional array of points to make a vector U of n² components). Then the relevant equation for each point U_i is:

(1) u_{x,x} + u_{y,y} = \frac{U_{i-1} + U_{i+1} +U_{i+n} + U_{i-n} - 4 U_{i}}{h^{2}} = k U_{i}

This is known to be true for this problem (see the link http://en.wikipedia.org/wiki/Discrete_Laplace_operator#Finite_Differences for the explanation)

and it results in a linear system of equations as follows:

(2) A U = kU + B

Where A is a sparse matrix which has five diagonals (1,1,-4,1,1) at the offsets ( -(n+1),-1,0,1,(n+1) ). Then B is another vector (of dimension n²) containing the boundary condition values - any point on the edge of the grid will contain a +1 term in equation (1) instead of the value of the point beside it (since U = 1 for all boundary points).

The system (2) should be solved by a linear solver (eg. conjugate gradient method or something similar) computationally. The result can be plotted as a surface and I believe all u(x,y) values should be between 0 and 1. (i.e. the surface should 'dip' down at all points within the boundaries.)

The Attempt at a Solution

[/B]
I understand the general theory behind writing the problem as a linear system of equations and that it should be solve-able as such. However when writing a program to solve this I'm unsure how to deal with

(a) The fact that the RHS of the linear system now depends on the solution U

(b) The boundary condition vector B

I'm used to solving linear systems of the form Ax=b where the right hand side is a constant, known vector. So my question is, how do I make the system (2) look like this? I will input a LHS and RHS but should all values of U be on the left, and which side does the vector B go on? Am I right in thinking that the RHS of the system must be constantly updated?

Any comments are appreciated, and I'm happy to clarify any other info about the problem. Thank you

Your linear equations are of the form
\frac{1}{h^2} (U_{i-1} + U_{i+1} +U_{i+n} + U_{i-n} - 4 U_{i}) = k U_i \\<br /> \;\;\;\;\;\;\text{or}\\<br /> \frac{1}{h^2}(U_{i-1} + U_{i+1} +U_{i+n} + U_{i-n}) - \left(\frac{4}{h^2} + k \right) U_i = 0 \\<br /> \;\;\;\;\;\; \text{or}\\<br /> U_{i-1} + U_{i+1} +U_{i+n} + U_{i-n} - (4 + k h^2) U_i = 0<br />

 

[beth92]

That makes sense - would it be correct then to say that: rather than necessarily having 0 on the right hand side for the above equation, I have B_{i} where B_{i} is:

-2 for the four points on the corner of the grid
-1 for points along the edge
0 for all other points

Thanks for the quick reply

 

[Ray Vickson]

That makes sense - would it be correct then to say that: rather than necessarily having 0 on the right hand side for the above equation, I have B_{i} where B_{i} is:

-2 for the four points on the corner of the grid
-1 for points along the edge
0 for all other points

Thanks for the quick reply

It is important to note that the difference-equation holds only for points not along an edge (because the pde holds in the interior of $\Omega$). So, in your 1-dimensional representation of the 2-dimensional region, let $S = \{1,2, \ldots N \}$ be the set of all possible $i$ values and $E \subset S$ be the points on an "edge". Then we ought to have
U_{i-1} + U_{i+1} +U_{i+n} + U_{i-n} - (4 + k h^2) U_i = 0, i \not \in E, \\<br /> U_i = 1, i \in E .
Of course, when $i \not \in E$, some of the $U_{i \pm 1}$ or $U_{i \pm n}$ could be in $E$, so you would substitute $1$ for them in that case. Basically, you have $N - |E|$ equations in $N - |E|$ unknowns (where $|E|$ denotes the cardinality of the set $E$).