1
$\begingroup$

I want to solve the following system via Matlab $\Omega=(0,1)^2$

$$\Delta y=\frac{1}{\alpha} p$$ $$ -\Delta p= y -1 $$ $$p|_{\partial \Omega}=0,~y|_{\partial \Omega}=0$$

using

    n=50;
    alpha=1;
    A=gallery('poisson',n);

    K=[A -1/alpha*speye(n^2);
       speye(n^2) A];
    b=[zeros(n^2,1) -ones(n^2,1)];
    solution=K\b;
    %embedding with boundaries...

To my surprise this leads to (I think) a correct answer, however I usually have this scaling factor $\frac{1}{h^2}=(n+1)^2$ in it. Taking $-\Delta y =1$ with zero boundaries as an example, this leads to $$\frac{1}{h^2}Ay=1_{n\times1}$$

Is my code incorrect, or does this factor vanish somehow?

$\endgroup$
  • 2
    $\begingroup$ Welcome to SciComp.SE! Are you sure the answer is right? I'd start with a manufactured solution (i.e., instead of $1$, pick a "nice" optimal state $y$ that satisfies the boundary conditions -- e.g., a product of sines and cosines -- compute (analytically!) the Laplacian to get the optimal adjoint state $p$, and then compute (analytically) the Laplacian and subtract $y$ to get the target $y_d$ (in place of $1$ here)); then you can see if the magnitude is correct. $\endgroup$ – Christian Clason Jan 16 '18 at 17:50
  • $\begingroup$ Thank you for your comment. I actually did this with $f(x,y)=sin(\pi x)sin(\pi y)$, and $y_d=f(x,y)(1-4 \alpha \pi^4 $ with the result $\|f-solution_h \|_\infty \le 10^{-5}$. For the proposed problem above, the outcome looks at least reasonable (correct is a bit strong). $\endgroup$ – user160069 Jan 16 '18 at 18:41
  • $\begingroup$ And does the error go to zero with the expected rate as $h\to 0$? Does it change when you change $\alpha$ (a lot)? $\endgroup$ – Christian Clason Jan 16 '18 at 18:55
  • $\begingroup$ choosing $\alpha=10^{-i}, i=2,3.4...150$ the error does get lower with increasing $i$, however for increasing $n$ it kind of stucks at this value. $\endgroup$ – user160069 Jan 16 '18 at 19:02
  • $\begingroup$ Then there is indeed a problem. What happens if you try with the correct scaling? $\endgroup$ – Christian Clason Jan 16 '18 at 19:16

Your Answer

By clicking "Post Your Answer", you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.