- #1
beth92
- 16
- 0
Homework Statement
I need to (computationally) solve the following linear elliptic problem for the function u(x,y):
[itex] \Delta u(x,y) = u_{x,x} + u_{y,y} = k u(x,y) [/itex]
on the domain
[itex] \Omega = [0,1]\times[0,1] [/itex] 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 n2 components). Then the relevant equation for each point Ui is:
(1) [itex] 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} [/itex]
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) [itex] A U = kU + B [/itex]
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 n2) 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 [itex]Ax=b[/itex] 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