Hi,

I am trying to solve the following 2D Poisson equation using the mglin function from numerical recipe in C (2nd Edition).

(d^2/dx^2 + d^2/dy^2)u = f
where f = -2*x*(1-x) - 2*y*(1-y)
& 0 < x,y < 2

Boundary conditions:
u = 0 for x=0, x=2, y=0, y=2

Following is my code:
-----------------------------------------------------------------------
#include<stdio.h>
#include<math.h>
#include "nrutil.c"
#include "nrutil.h"
#include "mglin.c"
#include "relax.c"
#include "slvsml.c"
#include "rstrct.c"
#include "copy.c"
#include "interp.c"
#include "fill0.c"
#include "resid.c"
#include "addint.c"

main()
{

int m = 1025;
double h = 2.0/(double)(m-1);

double **u;
u = dmatrix(1, m, 1, m);

int maxcyc = 2;

int i,j;
double x,y;


/*............Initializing u[i][j] = rhs[i][j].............*/

for(i=1;i<=m;i++)
for(j=1;j<=m;j++)
{
x = h*(i-1);
y = h*(j-1);

u[i][j] = -2*x*(1-x) -2*y*(1-y);

}

/************************************************** **********/



/*............Calling the mglin routine to solve the 2D PDE.....*/


mglin(u, m, maxcyc);

/************************************************** *****************/
}
----------------------------------------------------------------------

Changes made in slvsml.c, resid.c,relax.c:

1. slvsml.c :

In line 2: double h = 0.5 is replaced by h = 1.0

2. resid.c:

In line 3: double h = 1.0/(n-1) is replaced by 2.0/(n-1)

3. relax.c:

In line 3: double h = 1.0/(n-1) is replaced by 2.0/(n-1)

--------------------------------END--------------------------------------------

However, this does not give me the result which is x*(2-x)*y*(2-y)

Do you see where my mistake is?
Please let me know.

Thanks,
Arunima.

Dear Arunima,

The equation solved in NR has a u^2 term, which you don't have. So slvsml has a completely different solution, which you can easily work out.

Saul Teukolsky

Dear Sir,

I am using Numerical Recipe in C(2nd edition).
On pg.878, the mglin routine is given, which solves model problem 19.0.6
which is

(d^2/dx^2 + d^2/dy^2)u = rho

which is a linear problem. It doesnot contain u^2.

mgfas (pg.885) solves eq.19.6.44, which contains a u^2 term, and nonlinear.

But my equation is a linear one akin to 19.0.6, with the R.H.S a function of both the independent variables x & y.

Please let me know if I am mistaken. I am still confused.

Regards & thanks,
Arunima.

Sir,
I got it. The code is working fine with homogeneous dirichlet's b.c. I'll try my hands now on inhomogemeous/newman's b.c.
Thanks.

Arunima.

Hi,

I'm glad you got it to work. Sorry about the comment about slvsml - I misread your first post. You were correct that it did not need to be changed.

Saul Teukolsky

Hi,

I am trying to solve the following problem using mglin.

(d^2/dx^2 + d^2/dy^2)u = 2 (0 < = x,y < =1)
with boundary conditions as follows:

1. At x=0, u=0
2. At x=1, u=1
3. At y=0, du/dy = 0
4. At y=1, du/dy = 0

Following is my code:
-----------------------------------
#include<stdio.h>
#include<math.h>
#include "nrutil.c"
#include "nrutil.h"
#include "mglin.c"
#include "relax.c"
#include "slvsml.c"
#include "rstrct.c"
#include "copy.c"
#include "interp.c"
#include "fill0.c"
#include "resid.c"
#include "addint.c"

main()
{

FILE *fp;

int m = 1025;
double h = 1.0/(m-1);

double **u;
u = dmatrix(1, m, 1, m);

int maxcyc = 2;

int i,j;
double x,y,z,answer;


/*.....................................Initializing u[i][j] = rhs[i][j]..........................................*/

for(i=1;i<=m;i++)
for(j=1;j<=m;j++)
{
u[i][j] = 2.0;

}



i = m-1;
for(j=1;j<=m;j++)
{

x = (i-1)*h;
y = (j-1)*h;
u[i][j] = - (1.0)/(h*h) + 2.0;

}






/************************************************** ************************************************** *************/



/*...... Calling the mglin routine to solve the 2D PDE.....................................*/

mglin(u, m, maxcyc);

/************************************************** ************************************************** ********************/
}
--------------------------------------------------------------------------------------------------------
No changes have been made in slvsml.c, resid.c, and relax.c

This is not giving the expected solution which is u = x^2.

Can you please figure out my mistake?
Also, I'd like to confirm if the mglin code given in NR in C assumes homogeneous boundary conditions for both Dirichlet and Neumann's boundary conditions.

Regards & thanks,
Arunima.

Hi,

Yes, the code assumes homogeneous Dirichlet boundary conditions. To treat inhomogeneous boundary conditions, use the method described to change the right-hand side. In the second edition, this discussion follows equation (19.4.15) of the text. (20.4.15 in the third edition.)

Saul Teukolsky