# scaling in discretized PDE system

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?

• 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. Jan 16 '18 at 17:50
• 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). Jan 16 '18 at 18:41
• And does the error go to zero with the expected rate as $h\to 0$? Does it change when you change $\alpha$ (a lot)? Jan 16 '18 at 18:55
• 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. Jan 16 '18 at 19:02
• Then there is indeed a problem. What happens if you try with the correct scaling? Jan 16 '18 at 19:16