slvsm2 in mgfas()
I am trying to modify slvsm2() in the mgfas routine.
I was confused when I checked the detail of slvsm2().
void slvsm2(double **u, double **rhs)
{
.......
fact=2.0/(h*h);
disc=sqrt(fact*fact+rhs[2][2]);
u[2][2] = -rhs[2][2]/(fact+disc);
}
I solved the second-degree equation and got u[2][2] = fact+disc,
not -rhs[2][2]/(fact+disc).
Why u[2][2] = -rhs[2][2]/(fact+disc) is the right answer??
Somebody have an idea ??
:confused:
Saul Teukolsky
09-22-2004, 06:13 PM
Multiply the numerator and the denominator by fact-disc and you'll get the form you got (you typed a + instead of a -). The reason the form in NR is better is that it avoids the cancellation error if rhs is small. This point is discussed further in section 5.6.
smagorinsky
06-09-2008, 11:27 PM
Could there be any strong reason for choosing "fact - disc" instead of its "+" counterpart as the answer for the model problem (19.6.44) and the discretized form (19.6.45)?
Saul Teukolsky
06-10-2008, 03:30 PM
As rho tends to 0, the solution must be u = 0. That determines which solution to use.
Saul Teukolsky
smagorinsky
06-10-2008, 09:15 PM
Thank you very much for the reply, Professor Teukolsky.
As rho tends to 0, the solution must be u = 0. That determines which solution to use.
It is actually the confusing point for me.
The first thing:
You do not mean "rho" as the R.H.S. of the equation (19.6.44), do you? The "rho" in (19.6.44) is prescribed on each grid and will not be changed during the solution, from my understanding. Am I wrong somewhere?
The second thing:
Unlike the standart FMG, in FAS algorithm "u" on the coarse grids is not the corrector itself (comparing (19.6.11) and (19.6.29)) and I cannot see the reason for "u" tending to 0.
Would you kindly clarify these issues?
Saul Teukolsky
06-11-2008, 02:15 PM
*If* rho were zero, the solution *would be* u=0. This is an analytic statement about the original problem, and has nothing to do directly with the variables in the code. You have to write the code to solve the original problem, so by imagining this particular case you can see which root of the quadratic is the valid one.
Saul Teukolsky
smagorinsky
06-11-2008, 08:02 PM
I see. Thank you so much.