Hi,

I'm an undergraduate student from the Urals State University and I was given a task involving solving non-linear elliptic PDE. So I've taken "Numerical Recipes" (my favorite book) and typed mgfas routine from it. Than I've written a simple main function, ehich looks as follows:

main() {
int i,j;
int n = 9, maxcyc = 2;
double **u = dmatrix(1, n, 1, n);
double **rho = dmatrix(1, n, 1, n);
double h = 1.0/(n-1), h2i = 1.0/(h*h);

// Right-hand side
for (i = 1; i <= n; i++)
for (j = 1; j <= n; j++)
rho[i][j] = 1.; // rho(x,y)=1.=const

copy(u,rho,n);
mgfas(u, n, maxcyc);

// displaying solution
for (i = 2; i < n; i++) {
for (j = 2; j < n; j++)
printf("%f\t", u[i][j]);
printf("\n");
}
printf("\n");

//displaying defect for eq. (19.6.45)
for (i = 2; i < n; i++) {
for (j = 2; j < n; j++)
printf("%f\t",
(h2i*(u[i+1][j]+u[i-1][j]+u[i][j+1]+u[i][j-1]-4.0*u[i][j])
+u[i][j]*u[i][j]-rho[i][j]));
printf("\n");
}
}

The output of a program was urprising for me (pay attention to defect values)

-0.017743 -0.027689 -0.032833 -0.034437 -0.032833 -0.027689 -0.017743
-0.027689 -0.044567 -0.053623 -0.056478 -0.053623 -0.044567 -0.027689
-0.032833 -0.053623 -0.065034 -0.068653 -0.065034 -0.053623 -0.032833
-0.034437 -0.056478 -0.068653 -0.072514 -0.068653 -0.056478 -0.034437
-0.032833 -0.053623 -0.065034 -0.068653 -0.065034 -0.053623 -0.032833
-0.027689 -0.044567 -0.053623 -0.056478 -0.053623 -0.044567 -0.027689
-0.017743 -0.027689 -0.032833 -0.034437 -0.032833 -0.027689 -0.017743

-0.001712 0.000000 -0.001711 0.000000 -0.001711 0.000000 -0.001712
0.000000 0.003168 0.000000 0.000000 0.000000 0.003168 0.000000
-0.001711 0.000000 0.001711 0.000000 0.001711 0.000000 -0.001711
0.000000 0.000000 0.000000 -0.006337 0.000000 0.000000 0.000000
-0.001711 0.000000 0.001711 0.000000 0.001711 0.000000 -0.001711
0.000000 0.003168 0.000000 0.000000 0.000000 0.003168 0.000000
-0.001712 0.000000 -0.001711 0.000000 -0.001711 0.000000 -0.001712


Some of equations of a system (19.6.45) are satisfied exactly, but the others are satisfied very roughly. Is it a known issue with mgfas(), or I've done something wrong?

Best regards,
Valentine Sinitsyn

I've solved the system (19.6.45) using Newton's method (again from NR book) and found an interesting thing. Although values for u[i][j] differed from the ones I've obtained by mgfas() by 10^(-6) or so, the defects for the system (19.6.45) were much more higher (as shown on my previous post). So I think implementing test case for u[] in a way I've done isn't the best thing. I call this "unstability" of eqs (19.6.45)...

Hi Val,

It seems to me that mgfas is working exactly as it should. It's designed to solve the equation to truncation error. The idea is that there's no point in solving the difference equations to a higher precision than the error in replacing the pde by a difference approximation. Instead of solving the equations with a different method as you did, you could have set e.g. ALPHA to 1e-6 and maxcyc to 6 to get the "exact" solution to the difference equations. But there's no point in doing this. (This is discussed in the Example Book under mgfas.)

Hope this helps make things clearer.

Thanks a lot for your reply. I've also found that setting NPRE and NPOST to larger values (I've used 10 and 10) solves the problem.