Differential Approximation with Boundary Conditions

[Thorra]

Hello! I have a nifty set of problems (or rather one problem, gradually building itself to be a great problem) that I like to collectively call "The final problem" as it is the last thing I need before I can take the exam in Numerical Methods.Information

There is given a Laplace equation $\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = 0$ in a quadratic-type area $R = {(x,y)|0<x<1;0<y<1}$ with Neuman and Robin type boundary conditions $u(x,0)=0$ ; $\frac{\partial u(x,1)}{\partial y}+\frac{2u(x,1)}{(1+x^2)+1} = \frac{1}{(1+x^2)+1}$ ; $u(0,y)=\frac{y}{1+y^2}$ ; $\frac{\partial u(1,y)}{\partial x} = - \frac{4y}{(4+y^2)^2}$ .
Problem's analytical solution is $u(x,y)=\frac{y}{1+x^2)+y^2}$.

Tasks to do:

1. Approximate (with 2nd approximation order) the equation on the lattice/model with an even step ($h$ and $g$) in each of the directions $x$ and $y$, splitting each of the intervals $0<x<1$; $0<y<1$ in according parts $N$ and $M$.
2. Approximate (with 2nd approximation order) the given boundary conditions. For Neuman and Robin boundary conditions use fictive "nodal points" outside the area.

Okay... This is actually just the tip of the iceberg (apparently later I'll have to do programming... haha) but I already have trouble with this.

So I guess the 1st one I've kinda done before. I think it's something like:
$$\frac{u(x+h,y)-2u(x,y)+u(x-h,y)}{h^2}+\frac{u(x,y+g)-2u(x,y)+u(x,y-g)}{g^2}=0$$
where I $h$ and $g$ are divided into N and M parts $h=\frac{1}{N}$ and $g=\frac{1}{M}$ (1 being the interval length).
Errors are the old $\psi=\mathcal O(h^2+g^2)$.Now the 2nd one is a little trickier to me.
I counted 4 boundary conditions. Started with this one:
$$\frac{\partial u(1,y)}{\partial x}=-\frac{4y}{(4+y^2)^2}$$
I did this off a very similar example, but still not sure about it and not comprehending everything.
$$u(1+h,y)=u(1,y) + h\frac{\partial u(1,y)}{\partial x} + \frac{h^2}{2}\frac{\partial^2 u(1,y)}{\partial x^2} + O(h^3)$$
So that's just the ol' Taylor formula with the twist of taking the end point of x and stepping out of the model (this is what I naively interpret as the "fictive nodal point". It says to only use them on "Neuman" and "Robin" boundaries, but I don't think there even are any other boundary conditions).
Next we do the approximation of a first order derivative:
$$\frac{u(1+h,y)-u(1,y)}{h}=\frac{\partial u(1,y)}{\partial x}+\frac{h}{2}\frac{\partial^2 u(1,y)}{\partial x^2}+\mathcal O(h^2)$$
And now a really confusing moment for me:
$$\frac{u(1+h,y)-u(1,y)}{h}=-\frac{4y}{(4+y^2)^2}-\frac{h}{2}\frac{\partial^2 u(x,y)}{\partial y^2}+\mathcal O(h^2)$$
It's not the small substitution at the first part. It's the 2nd part where the u(1,y) magically turns into u(x,y) and x^2 turns into y^2. I just don't understand why would that be. It basically means that $$\frac{\partial^2 u(1,y)}{\partial x^2}=-\frac{\partial^2 u(x,y)}{\partial y^2}$$Thanks in advance.

 

[Thorra]

Okay... I've done quite some research and I can only conclude that:
$u(x,0)=0$, $u(0,y)=\frac{y}{1+y^2}$ - Dirichlet boundary conditions (why imply that there are only Neumann and Robin BCs? Wtf..)
$\frac{\partial u(x,1)}{\partial y}+\frac{2u(x,1)}{(1+x^2)+1} = \frac{1}{(1+x^2)+1}$ - Robin BC
$\frac{\partial u(1,y)}{\partial x} = - \frac{4y}{(4+y^2)^2}$ - Neumann BC.

Can anyone confirm this?

 

[I like Serena]

Hello! I have a nifty set of problems (or rather one problem, gradually building itself to be a great problem) that I like to collectively call "The final problem" as it is the last thing I need before I can take the exam in Numerical Methods.Information

There is given a Laplace equation $\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = 0$ in a quadratic-type area $R = {(x,y)|0<x<1;0<y<1}$ with Neuman and Robin type boundary conditions $u(x,0)=0$ ; $\frac{\partial u(x,1)}{\partial y}+\frac{2u(x,1)}{(1+x^2)+1} = \frac{1}{(1+x^2)+1}$ ; $u(0,y)=\frac{y}{1+y^2}$ ; $\frac{\partial u(1,y)}{\partial x} = - \frac{4y}{(4+y^2)^2}$ .
Problem's analytical solution is $u(x,y)=\frac{y}{1+x^2)+y^2}$.

There's a couple of apparent inconsistencies.
For starters, there's a parenthesis missing.

And if I substitute the analytical solution, the 2nd boundary condition does not match.
Perhaps a square was dropped somewhere?
It does look strange to write $(1+x^2)+y^2$.

Tasks to do:

1. Approximate (with 2nd approximation order) the equation on the lattice/model with an even step ($h$ and $g$) in each of the directions $x$ and $y$, splitting each of the intervals $0<x<1$; $0<y<1$ in according parts $N$ and $M$.

Okay... This is actually just the tip of the iceberg (apparently later I'll have to do programming... haha) but I already have trouble with this.

So I guess the 1st one I've kinda done before. I think it's something like:
$$\frac{u(x+h,y)-2u(x,y)+u(x-h,y)}{h^2}+\frac{u(x,y+g)-2u(x,y)+u(x,y-g)}{g^2}=0$$
where I $h$ and $g$ are divided into N and M parts $h=\frac{1}{N}$ and $g=\frac{1}{M}$ (1 being the interval length).
Errors are the old $\psi=\mathcal O(h^2+g^2)$.

Looks good. (Nod)

Actually, I consider the programming the easy part. (Angel)

2. Approximate (with 2nd approximation order) the given boundary conditions. For Neuman and Robin boundary conditions use fictive "nodal points" outside the area.

Now the 2nd one is a little trickier to me.
I counted 4 boundary conditions. Started with this one:
$$\frac{\partial u(1,y)}{\partial x}=-\frac{4y}{(4+y^2)^2}$$
I did this off a very similar example, but still not sure about it and not comprehending everything.
$$u(1+h,y)=u(1,y) + h\frac{\partial u(1,y)}{\partial x} + \frac{h^2}{2}\frac{\partial^2 u(1,y)}{\partial x^2} + O(h^3)$$
So that's just the ol' Taylor formula with the twist of taking the end point of x and stepping out of the model (this is what I naively interpret as the "fictive nodal point". It says to only use them on "Neuman" and "Robin" boundaries, but I don't think there even are any other boundary conditions).
Next we do the approximation of a first order derivative:
$$\frac{u(1+h,y)-u(1,y)}{h}=\frac{\partial u(1,y)}{\partial x}+\frac{h}{2}\frac{\partial^2 u(1,y)}{\partial x^2}+\mathcal O(h^2)$$
And now a really confusing moment for me:
$$\frac{u(1+h,y)-u(1,y)}{h}=-\frac{4y}{(4+y^2)^2}-\frac{h}{2}\frac{\partial^2 u(x,y)}{\partial y^2}+\mathcal O(h^2)$$
It's not the small substitution at the first part. It's the 2nd part where the u(1,y) magically turns into u(x,y) and x^2 turns into y^2. I just don't understand why would that be. It basically means that $$\frac{\partial^2 u(1,y)}{\partial x^2}=-\frac{\partial^2 u(x,y)}{\partial y^2}$$

That looks like a typo.

From the Laplace equation, we have:
\begin{array}{}
\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} &=& 0 & \qquad (1) \\
\frac{\partial^2 u}{\partial x^2} &=& -\frac{\partial^2 u}{\partial y^2} \\
\frac{\partial^2 u(1,y)}{\partial x^2} &=& -\frac{\partial^2 u(1,y)}{\partial y^2}
\end{array}

Okay... I've done quite some research and I can only conclude that:
$u(x,0)=0$, $u(0,y)=\frac{y}{1+y^2}$ - Dirichlet boundary conditions (why imply that there are only Neumann and Robin BCs? Wtf..)
$\frac{\partial u(x,1)}{\partial y}+\frac{2u(x,1)}{(1+x^2)+1} = \frac{1}{(1+x^2)+1}$ - Robin BC
$\frac{\partial u(1,y)}{\partial x} = - \frac{4y}{(4+y^2)^2}$ - Neumann BC.

Can anyone confirm this?

I have just looked them up myself and I draw the same conclusions.
Perhaps other people can make mistakes as well. ;)
Luckily, it doesn't matter for the problem at hand.

 

[Thorra]

There's a couple of apparent inconsistencies.
For starters, there's a parenthesis missing.
More importantly, it makes no sense to write $(1+x^2)+1$.
Can it be that a square was dropped or something like that?

Nope, I rechecked and it is written that way. I know, weird, but maybe it's some implication from where it came from or whatever.

Actually, I consider the programming the easy part. (Angel)

Ha, I wish I could say that. I have 8 points on there and I want at least 4 done by monday so I can hope for a passing grade. Still stuck at 2!

That looks like a typo.

From the Laplace equation, we have:
\begin{array}{}
\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} &=& 0 & \qquad (1) \\
\frac{\partial^2 u}{\partial x^2} &=& -\frac{\partial^2 u}{\partial y^2} \\
\frac{\partial^2 u(1,y)}{\partial x^2} &=& -\frac{\partial^2 u(1,y)}{\partial y^2}
\end{array}

Ooooh! Well that solves one riddle.

I have just looked them up myself and I draw the same conclusions.
Perhaps other people can make mistakes as well. ;)
Luckily, it doesn't matter for the problem at hand.

Okay, but I must ask, not sure if you know. I read in a book that with the Robin and Neuman BCs you pretty much have to introduce the fictive nodes (to save the matrix structure and 2nd order of approximation), so could I have it right altogether? I did make one mistake for sure at the beginning though I think.

rather than $u(1+h,y)=u(1,y) +...$ I should probably choose $u(1-h,y)$ to stay within the model.

So yeah, that's it more or less?

Since you consider programing so easy, how what do you think bout the task of:
4. Solve the acquired Linear Algebraic Equation System (LAES), using any of Gauss exclusion technique algorithms and sequentially taking $h=g$ and $N=M=4;6;8;10;12;14;16$ .

Admittely I went straight over a task to write out the N=M=4 matrix out but I'll try to figure that one out myself.

I don't know any algorithms and it looks very programmy, all I can tell is that it would have an equally tall-long matrix and it would be 16x16, 36x36, etc. if I'm even correct.

Luckily, it doesn't matter for the problem at hand.

Why does it not matter?? I don't comprehend.

 

[I like Serena]

Nope, I rechecked and it is written that way. I know, weird, but maybe it's some implication from where it came from or whatever.

The 2nd boundary condition does not match with the analytic solution...

Okay, but I must ask, not sure if you know. I read in a book that with the Robin and Neuman BCs you pretty much have to introduce the fictive nodes (to save the matrix structure and 2nd order of approximation), so could I have it right altogether? I did make one mistake for sure at the beginning though I think.

Your 2nd equation depends on a 2nd derivative.
It seems to me that should be eliminated.

Or rather, if it were up to me, I'd use the approximation:
$$f'(x) = \frac{f(x+h) - f(x-h)}{2h} + \mathcal O(h^2)$$

rather than $u(1+h,y)=u(1,y) +...$ I should probably choose $u(1-h,y)$ to stay within the model.

Is that important?
You might also introduce fictive nodes.
Unless there is reason to believe that the function breaks down somewhere outside the given domain.

Since you consider programing so easy, how what do you think bout the task of:
4. Solve the acquired Linear Algebraic Equation System (LAES), using any of Gauss exclusion technique algorithms and sequentially taking $h=g$ and $N=M=4;6;8;10;12;14;16$.

Admittely I went straight over a task to write out the N=M=4 matrix out but I'll try to figure that one out myself.

I don't know any algorithms and it looks very programmy, all I can tell is that it would have an equally tall-long matrix and it would be 16x16, 36x36, etc. if I'm even correct.

With N=M=4 you would get 5x5 matrices (including the boundaries).
Or up to 7x7 matrices if you include fictive nodes.

Why does it not matter?? I don't comprehend.

What's in a name?
So one of the conditions is not actually a Neuman or Robin condition.
For the problem it's only the equation of the boundary that matters, not how it is called.

It does matter though in whether you may want to introduce fictive nodes or not on the particular boundary.

 

[Thorra]

The 2nd boundary condition does not match with the analytic solution...

is that bad? With my level of underswtanding I think only the simple boundary conditions match the analytic solution [u(x,0) and u(0,y)]

Your 2nd equation depends on a 2nd derivative.
It seems to me that should be eliminated.

Or rather, if it were up to me, I'd use the approximation:
$$f'(x) = \frac{f(x+h) - f(x-h)}{2h} + \mathcal O(h^2)$$

Which 2nd equation do you mean? Do you mean like this?
$$-\frac{4y}{4+y^2)^2}=\frac{u(1+h) - u(1-h)}{2h} + \mathcal O(h^2)$$
sorry if I'm shooting really wrong, my brain is not what it was in the morning and I am really in a hurry here so I can pass the subject.

Is that important?
You might also introduce fictive nodes.
Unless there is reason to believe that the function breaks down somewhere outside the given domain.

Yeah that's what I'm asked to do and I don't know how to do. I read all kinds of theory and looked at pretty graphs that described my situation, but nothing about actual approximation with them.

With N=M=4 you would get 5x5 matrices (including the boundaries).
Or up to 7x7 matrices if you include fictive nodes.

Really? I was lead to believe that it'd be a 16x16 matrix consisting of 16 4x4 submatrix blocks. Without fictive nodes.

What's in a name?
So one of the conditions is not actually a Neuman or Robin condition.
For the problem it's only the equation of the boundary that matters, not how it is called.

I thought it changes something in the way you use/approximate them otherwise why bother to sepcify them.

It does matter though in whether you may want to introduce fictive nodes or not on the particular boundary.

I do! But how? I was told that approximating with a step out of the model isn't the way to do it which was my intuitive thought (haha, yeah, using intuition in this subject of hell is not a thing to do).

Or maybe I use something like $$\frac{u(1+2h,y)-u(1,y)}{h}$$ ?EDIT: I guess that didn't do much... Ended up with the same thing.
$$\omega^j_{N+2}-\omega^j_N = -\frac{4y}{(4+y^2)^2}-\frac{h}{2}\frac{\omega^{j+1}_N - 2\omega^j_N+\omega^{j-1}_N}{g^2}$$

 

[I like Serena]

is that bad? With my level of underswtanding I think only the simple boundary conditions match the analytic solution [u(x,0) and u(0,y)]

The analytic solution is only a solution if it satisfies the Laplace equation and all boundary conditions. Since it does not satisfy the second boundary condition, it is not a solution.

In other words, I think there is a typo somewhere.
It would be good to find it and eliminate it.

Which 2nd equation do you mean? Do you mean like this?
$$-\frac{4y}{4+y^2)^2}=\frac{u(1+h) - u(1-h)}{2h} + \mathcal O(h^2)$$
sorry if I'm shooting really wrong, my brain is not what it was in the morning and I am really in a hurry here so I can pass the subject.

Yes (after you fix the typos and ambiguities).

Yeah that's what I'm asked to do and I don't know how to do. I read all kinds of theory and looked at pretty graphs that described my situation, but nothing about actual approximation with them.

Really? I was lead to believe that it'd be a 16x16 matrix consisting of 16 4x4 submatrix blocks. Without fictive nodes.

Well, perhaps if you write it out they will be.
But then you will first have to write out at least 1 example.

I thought it changes something in the way you use/approximate them otherwise why bother to sepcify them.

Specify what? The names or the boundary conditions?
Giving names to stuff helps to sort it out.
At least now you know what Direchlet, Neumann, and Robin conditions are.
And actually they each have different requirements for how to deal with them. For instance by using either fictive nodes or not.

I do! But how? I was told that approximating with a step out of the model isn't the way to do it which was my intuitive thought (haha, yeah, using intuition in this subject of hell is not a thing to do).

Intuition is good! Of course it has to be backed up by numbers.

Or maybe I use something like $$\frac{u(1+2h,y)-u(1,y)}{h}$$ ?EDIT: I guess that didn't do much... Ended up with the same thing.
$$\omega^j_{N+2}-\omega^j_N = -\frac{4y}{(4+y^2)^2}-\frac{h}{2}\frac{\omega^{j+1}_N - 2\omega^j_N+\omega^{j-1}_N}{g^2}$$

Huh? I don't understand.

It seems to me you need a matrix with values for $u(x + jh, y + ig)$ where $j=-1,0,..., N+1$ and $i=-1,0,...,M+1$.
This is the "solution" to your system.
It has to satisfy (the approximation of) the Laplace equation and also the boundary conditions.

 

[Thorra]

The analytic solution is only a solution if it satisfies the Laplace equation and all boundary conditions. Since it does not satisfy the second boundary condition, it is not a solution.

In other words, I think there is a typo somewhere.
It would be good to find it and eliminate it.

Hmm.. maybe adding the quadratic expression after the parenthesis.
$\frac{\partial u(x,1)}{\partial y}+\frac{2u(x,1)}{(1+x)^2+1} = \frac{1}{(1+x)^2+1}$
That would explain why have separate brackets for it. :D

Yes (after you fix the typos and ambiguities).

Like this, then?
$$-\frac{4y}{(4+y^2)^2}=\frac{u(1+h,y) - u(1-h,y)}{2h} + \mathcal O(h^2)$$

Well, perhaps if you write it out they will be.
But then you will first have to write out at least 1 example.

Yeah and that's probably why I need to find all the right expressions first...
So I'm trying to wrap my head around this.
1. I need a lot of values for the matrix. For each point. I got the ones specified with the bottom BCs known... and through something I cannot yet accomplish I can find out the upper boundaries in a linear, useful way so that I could put that stuff in the matrixes.
2. I think I got this situation visually:

I think that because it also has a dirichlet condition on the bottom x and y limits and a neuman for upper x limit and robin for upper y limit.
So I need to find all the white points, with the help of fictious nodes and (probably) black nodes. All points from j=1 to j=N=4 and k=1 to k=M=4.
3. So how do I do that? The book mentions this equation: $\frac{u_{j-1,k}-2u_{j,k}+u{j+1,k}}{h^2}+\mathcal O(h^2)+\frac{u_{j-1,k}-2u_{j,k}+u{j+1,k}}{h^2}+\mathcal O(h^2)=f_{j,k}$ that I multiply by $(-h^2)$ and get a nice function $-\omega_{j-1,k}-\omega_{j+1,k}-\omega_{j,k-1}-\omega{j,k+1}+4\omega_{j,k}=-h^2 f_{j,k}$. Sorry for the long typing, but I don't know what other direction to follow. And I can't even follow this! I don't have equivalent steps. Or maybe I do, but it's not specified. So what do I do?
Later it mentions how I can substitute Dirichlet boundary territory (black nodes) $\omega_{1,0}=g_{1,0}$ and $\omega_{0,1}=g_{0,1}$ and bring them to the right side of the equation (the one that is not possible in that form for me). Like this: $4\omega_{1,1}-\omega_{2,1}-\omega_{1,2}=-h^2f_{1,1}+g_{1,0}+g_{0,1}$.

I'm sorry, but I don't know what I'm doing anymore. The more I go into this stuff the less of an idea I have what I need to achieve and how to achieve anything.

There is also this thing of the Robin BC (the most complicated one) but I don't expect you to make much sense of it: $$(-1-\frac{h}{2}p_i)\omega_{i-1}+(2+h^2q_i)\omega_i+(-1+\frac{h}{2}p_i)\omega_{i+1}=-h^2r_i$$
where you make i=0 and get in a fictional value $\omega_f$ in the game and do more things with that. Maybe that's the way to go? It's only for x-axis though, it seems, so I hope that's not a big bother.

Specify what? The names or the boundary conditions?
Giving names to stuff helps to sort it out.
At least now you know what Direchlet, Neumann, and Robin conditions are.
And actually they each have different requirements for how to deal with them. For instance by using either fictive nodes or not.

Yeah... that's what I mean. I want to use those fictive values.
For Neumann's:$$\frac{\omega_f-\omega_N-1}{2h}=-\frac{4y}{(4+y^2)^2}$$ $$\omega_f=\omega_N-1-\frac{4y}{(4+y^2)^2}$$
and put that in... somewhere... and get a value at a right-side (of the drawing) boundary point!
Somewhere? Mm... maybe...the original equation?
$$\frac{\omega(x+h,y)-2\omega(x,y)+\omega(x-h,y)}{h^2}+\frac{\omega(x,y+g)-2\omega(x,y)+\omega(x,y-g)}{g^2}=0$$
$$\frac{\omega_f^j-2\omega_N^j+\omega_{N-1}^j}{h^2}$$
Haha.. notice how I don't know what to do with the other part. I hate this subject so much. Thank you for sticking with me through it, man.

Huh? I don't understand.

It seems to me you need a matrix with values for $u(x + jh, y + ig)$ where $j=-1,0,..., N+1$ and $i=-1,0,...,M+1$.
This is the "solution" to your system.
It has to satisfy (the approximation of) the Laplace equation and also the boundary conditions.

x+jh? Don't you mean x+ih? I hope to god you just innocently mixed that up. And why j=-1 and i=-1? I think the bottom and left side are restricted to the known dirichlet conditions where fictious nodes are not needed.

 

[I like Serena]

Hmm.. maybe adding the quadratic expression after the parenthesis.
$\frac{\partial u(x,1)}{\partial y}+\frac{2u(x,1)}{(1+x)^2+1} = \frac{1}{(1+x)^2+1}$
That would explain why have separate brackets for it. :D

I checked the 2nd boundary condition again and it turns out I made a mistake.
The analytic solution does satisfy it.
So there should only be an extra parenthesis in the analytic solution.

Like this, then?
$$-\frac{4y}{(4+y^2)^2}=\frac{u(1+h,y) - u(1-h,y)}{2h} + \mathcal O(h^2)$$

Yep. ;)

Yeah and that's probably why I need to find all the right expressions first...
So I'm trying to wrap my head around this.
1. I need a lot of values for the matrix. For each point. I got the ones specified with the bottom BCs known... and through something I cannot yet accomplish I can find out the upper boundaries in a linear, useful way so that I could put that stuff in the matrixes.
2. I think I got this situation visually:

I think that because it also has a dirichlet condition on the bottom x and y limits and a neuman for upper x limit and robin for upper y limit.
So I need to find all the white points, with the help of fictious nodes and (probably) black nodes. All points from j=1 to j=N=4 and k=1 to k=M=4.

Aha! That is a very nice picture that shows exactly how your setup should be! (Cool)
So for M=N=4 you might use a 5x5 matrix to contain all the $u$-values.

3. So how do I do that? The book mentions this equation: $\frac{u_{j-1,k}-2u_{j,k}+u{j+1,k}}{h^2}+\mathcal O(h^2)+\frac{u_{j-1,k}-2u_{j,k}+u{j+1,k}}{h^2}+\mathcal O(h^2)=f_{j,k}$ that I multiply by $(-h^2)$ and get a nice function $-\omega_{j-1,k}-\omega_{j+1,k}-\omega_{j,k-1}-\omega{j,k+1}+4\omega_{j,k}=-h^2 f_{j,k}$. Sorry for the long typing, but I don't know what other direction to follow. And I can't even follow this! I don't have equivalent steps. Or maybe I do, but it's not specified. So what do I do?

What's the problem?
It is the Laplace equation with the 2nd derivatives replaced with approximations based on the $u$-values. In your case $f_{j,k}=0$.
What is the reason to switch to $\omega$ symbols instead of $u$ symbols?

Later it mentions how I can substitute Dirichlet boundary territory (black nodes) $\omega_{1,0}=g_{1,0}$ and $\omega_{0,1}=g_{0,1}$ and bring them to the right side of the equation (the one that is not possible in that form for me). Like this: $4\omega_{1,1}-\omega_{2,1}-\omega_{1,2}=-h^2f_{1,1}+g_{1,0}+g_{0,1}$.
I'm sorry, but I don't know what I'm doing anymore. The more I go into this stuff the less of an idea I have what I need to achieve and how to achieve anything.

Same thing.
The Dirichlet boundary conditions are written using $u$-values from a matrix.
It's only confusing that suddenly $\omega$ is used instead of $u$.

There is also this thing of the Robin BC (the most complicated one) but I don't expect you to make much sense of it: $$(-1-\frac{h}{2}p_i)\omega_{i-1}+(2+h^2q_i)\omega_i+(-1+\frac{h}{2}p_i)\omega_{i+1}=-h^2r_i$$
where you make i=0 and get in a fictional value $\omega_f$ in the game and do more things with that. Maybe that's the way to go? It's only for x-axis though, it seems, so I hope that's not a big bother.

This looks like a generic way to write a Robin BC with discrete values and approximations.
You're supposed to write your own Robin BC in that form.

Yeah... that's what I mean. I want to use those fictive values.
For Neumann's:$$\frac{\omega_f-\omega_N-1}{2h}=-\frac{4y}{(4+y^2)^2}$$ $$\omega_f=\omega_N-1-\frac{4y}{(4+y^2)^2}$$
and put that in... somewhere... and get a value at a right-side (of the drawing) boundary point!
Somewhere? Mm... maybe...the original equation?
$$\frac{\omega(x+h,y)-2\omega(x,y)+\omega(x-h,y)}{h^2}+\frac{\omega(x,y+g)-2\omega(x,y)+\omega(x,y-g)}{g^2}=0$$
$$\frac{\omega_f^j-2\omega_N^j+\omega_{N-1}^j}{h^2}$$
Haha.. notice how I don't know what to do with the other part. I hate this subject so much. Thank you for sticking with me through it, man.

More of the same...
... but what's $f$ supposed to represent?

x+jh? Don't you mean x+ih? I hope to god you just innocently mixed that up.

Neh, not innocently. (Angel)

It's a matter of preference.
In a matrix rows are often identified by $i$ and the columns by $j$.
I like to have the x-coordinate horizontal and the y-coordinate vertical.
To obtain that, $i$ should be the index for the y-coordinate and $j$ the index for the x-coordinate.

It appears that in your book example the choices for indices and coordinates was made differently.

And why j=-1 and i=-1? I think the bottom and left side are restricted to the known dirichlet conditions where fictious nodes are not needed.

You're quite right.
And even for the Neumann and Robin BC's you have a choice.
You can either use an approximation that uses fictitious points, or you can use an approximation for the derivative that only uses points inside and on the boundary of the domain.

 

[Thorra]

I checked the 2nd boundary condition again and it turns out I made a mistake.
The analytic solution does satisfy it.
So there should only be an extra parenthesis in the analytic solution.

Like this you mean? $u(x,y)=\frac{y}{(1+x^2)+y^2}$.
That's just a minor typo if that's what you meant.

Aha! That is a very nice picture that shows exactly how your setup should be! (Cool)
So for M=N=4 you might use a 5x5 matrix to contain all the $u$-values.

So what should the matrix look like? Because I cannot for the life of me figure out how it works for real. The books are useless.
I can only think of this old thing (sorry, sticking with the old system of i's and j's. It's simpler for me for now)
$$\frac{1}{h^2}\omega^j_{i+1}+\frac{1}{g^2}\omega^{j+1}_i+\left(-\frac{2}{h^2}-\frac{2}{g^2}\right)+\frac{1}{h^2}\omega^j_{i-1}+\frac{1}{g^2}\omega^{j-1}_i=0$$
and I put in there i=0,j=0 and i=1,j=1 up to i,j=4 and I assign the appropriate $\omega$ coefficients to the appropriate matrix fields and get...
\[\left(\begin{matrix}\left(-\frac{2}{h^2}-\frac{2}{g^2}\right)&\frac{1}{g^2}&0?&0?&0?\\\frac{1}{h^2}&\left(-\frac{2}{h^2}-\frac{2}{g^2}\right)&\frac{1}{g^2}&0?&0?\\0&\frac{1}{h^2}&\left(-\frac{2}{h^2}-\frac{2}{g^2}\right)&\frac{1}{h^2}&0\\0&0&\frac{1}{h^2}&\left(-\frac{2}{h^2}-\frac{2}{g^2}\right)&\frac{1}{g^2}\\0&0&0&\frac{1}{h^2}&\left(-\frac{2}{h^2}-\frac{2}{g^2}\right)\end{matrix}\right)\]
I know it's completely wrong, but I cannot figure out anything else. Could you throw any idea or magic formulas that would fix everything? It may look like I'm not even trying, but the sad part is I am trying my best to dig through this, spending all my time on it.

What's the problem?
It is the Laplace equation with the 2nd derivatives replaced with approximations based on the $u$-values. In your case $f_{j,k}=0$.
What is the reason to switch to $\omega$ symbols instead of $u$ symbols?

I switch to $\omega$ because I drop the $\mathcal O(h^2)$ and instead use an approximate value of $u$.
And the problem is that I shouldn't just assume that my g=h (even tho it probably is), so I can't tell out the $\omega_i$ I'm searching for.

This looks like a generic way to write a Robin BC with discrete values and approximations.
You're supposed to write your own Robin BC in that form.More of the same...
... but what's $f$ supposed to represent?

Well this was me writing "my own" Neumann BC. I thought that's the way to do it. $f$ is supposed to mean fictious! For the fictious nodes.

 

[I like Serena]

Like this you mean? $u(x,y)=\frac{y}{(1+x^2)+y^2}$.
That's just a minor typo if that's what you meant.

Yes. It turns out it is just a minor typo.
It was bugging me though that there appeared to be serious typos in the problem statement.
It's a good thing there aren't. :)

So what should the matrix look like? Because I cannot for the life of me figure out how it works for real. The books are useless.
I can only think of this old thing (sorry, sticking with the old system of i's and j's. It's simpler for me for now)
$$\frac{1}{h^2}\omega^j_{i+1}+\frac{1}{g^2}\omega^{j+1}_i+\left(-\frac{2}{h^2}-\frac{2}{g^2}\right)+\frac{1}{h^2}\omega^j_{i-1}+\frac{1}{g^2}\omega^{j-1}_i=0$$
and I put in there i=0,j=0 and i=1,j=1 up to i,j=4 and I assign the appropriate $\omega$ coefficients to the appropriate matrix fields and get...
\[\left(\begin{matrix}\left(-\frac{2}{h^2}-\frac{2}{g^2}\right)&\frac{1}{g^2}&0?&0?&0?\\\frac{1}{h^2}&\left(-\frac{2}{h^2}-\frac{2}{g^2}\right)&\frac{1}{g^2}&0?&0?\\0&\frac{1}{h^2}&\left(-\frac{2}{h^2}-\frac{2}{g^2}\right)&\frac{1}{h^2}&0\\0&0&\frac{1}{h^2}&\left(-\frac{2}{h^2}-\frac{2}{g^2}\right)&\frac{1}{g^2}\\0&0&0&\frac{1}{h^2}&\left(-\frac{2}{h^2}-\frac{2}{g^2}\right)\end{matrix}\right)\]
I know it's completely wrong, but I cannot figure out anything else. Could you throw any idea or magic formulas that would fix everything? It may look like I'm not even trying, but the sad part is I am trying my best to dig through this, spending all my time on it.

Let's call the matrix with $u$ values $U$.
If you multiply $U$ on the right side with:
\[H = \left(\begin{matrix}\left(-\frac{2}{h^2}\right)&\frac{1}{h^2}&0&0&0\\
\frac{1}{h^2}&\left(-\frac{2}{h^2}\right)&\frac{1}{h^2}&0&0\\
0&\frac{1}{h^2}&\left(-\frac{2}{h^2}\right)&\frac{1}{h^2}&0\\
0&0&\frac{1}{h^2}&\left(-\frac{2}{h^2}\right)&\frac{1}{h^2}\\
0&0&0&\frac{1}{h^2}&\left(-\frac{2}{h^2}\right)\end{matrix}\right)\]
you get the matrix with the \(\displaystyle \frac{\partial^2 u(x,y)}{\partial x^2}\) values.

Similarly multiply $U$ on the left side with a matrix with the $g$ values (which we'll call $G$), and you will get the \(\displaystyle \frac{\partial^2 u(x,y)}{\partial y^2}\) values.

In other words, your Laplace equation is:
$$UH + GU = 0$$

Similarly, you can define the boundary conditions.

I switch to $\omega$ because I drop the $\mathcal O(h^2)$ and instead use an approximate value of $u$.
And the problem is that I shouldn't just assume that my g=h (even tho it probably is), so I can't tell out the $\omega_i$ I'm searching for.

Ah okay.
And indeed, you shouldn't assume g=h.

Well this was me writing "my own" Neumann BC. I thought that's the way to do it. $f$ is supposed to mean fictious! For the fictious nodes.

Hmm, as I see it, those fictitious nodes should have normal $i$ and $j$ indices.
Just indices that are "out-of-range". So with i=M+1 you would address a fictitious node.

 

[Thorra]

Let's call the matrix with $u$ values $U$.
If you multiply $U$ on the right side with:
\[H = \left(\begin{matrix}\left(-\frac{2}{h^2}\right)&\frac{1}{h^2}&0&0&0\\
\frac{1}{h^2}&\left(-\frac{2}{h^2}\right)&\frac{1}{h^2}&0&0\\
0&\frac{1}{h^2}&\left(-\frac{2}{h^2}\right)&\frac{1}{h^2}&0\\
0&0&\frac{1}{h^2}&\left(-\frac{2}{h^2}\right)&\frac{1}{h^2}\\
0&0&0&\frac{1}{h^2}&\left(-\frac{2}{h^2}\right)\end{matrix}\right)\]
you get the matrix with the \(\displaystyle \frac{\partial^2 u(x,y)}{\partial x^2}\) values.

Similarly multiply $U$ on the left side with a matrix with the $g$ values (which we'll call $G$), and you will get the \(\displaystyle \frac{\partial^2 u(x,y)}{\partial y^2}\) values.

In other words, your Laplace equation is:
$$UH + GU = 0$$

Similarly, you can define the boundary conditions.

So
$U\cdot(H+G)=$\begin{matrix}\left(-\frac{2}{h^2}-\frac{2}{g^2}\right)&\frac{1}{h^2}+\frac{1}{g^2}&0&0&0\\
\frac{1}{h^2}+\frac{1}{g^2}&\left(-\frac{2}{h^2}-\frac{2}{g^2}\right)&\frac{1}{h^2}+\frac{1}{g^2}&0&0\\
0&\frac{1}{h^2}+\frac{1}{g^2}&\left(-\frac{2}{h^2}-\frac{2}{g^2}\right)&\frac{1}{h^2}+\frac{1}{g^2}&0\\
0&0&\frac{1}{h^2}+\frac{1}{g^2}&\left(-\frac{2}{h^2}-\frac{2}{g^2}\right)&\frac{1}{h^2}+\frac{1}{g^2}\\
0&0&0&\frac{1}{h^2}+\frac{1}{g^2}&\left(-\frac{2}{h^2}-\frac{2}{g^2}\right)\end{matrix}
So it's that simple? I saw a big and complicated 16x16 matrix for this same equation.
To work in the boundary conditions though I'll have to figure out the boundary code and approximations and whatnot. Are boundary conditions the ones the complicate the matrxi from a 5x5 (or 4x4 if you don't count the black nodes or something) to 16x16? I don't know I'll think about it.

 

[I like Serena]

Sorry, you can't change the order of matrix multiplication.
In general:
$$GU + UH \ne UG + UH$$

 

[Thorra]

I have a question: what do I need to change in this process for it to have approximation with fictive nodes? Rather than non fictive nodes. (It's both the Neuman and Robin BCs approximated with I don't know what method.)

 

[I like Serena]

I have a question: what do I need to change in this process for it to have approximation with fictive nodes? Rather than non fictive nodes. (It's both the Neuman and Robin BCs approximated with I don't know what method.)

Sorry, but those are a couple too many formulas with too little context for me to give an answer.

 

[Thorra]

Hello again.

I got some help with this since I last posted, but I ran into a problem later on (when calculating all the points on the models when N=M=4,6,8,..,16) which suggested the whole job was done with the 1st level of approximation. So could the programming method (yes, there was programming) be responsible for that? (It was not done with iterations or anything, like some coursemates did it) Or could I have approximated the essential scheme wrong?

Same boundary condition
$$\frac{\partial u(1,y)}{\partial x}=-\frac{4y}{(4+y^2)^2}$$
approximated with central difference.
$$u(1-h,y)=u(1,y)-h\frac{\partial u(1,y)}{\partial x}+\frac{h^2}{2}\frac{\partial^2 u(1,y)}{\partial x^2}+\mathcal O(h^3)$$
$$u(1+h,y)=u(1,y)+h\frac{\partial u(1,y)}{\partial x}+\frac{h^2}{2}\frac{\partial^2 u(1,y)}{\partial x^2}+\mathcal O(h^3)$$
$$\frac{u(1+h,y)-u(1-h,y)}{2h}=\frac{\partial u(1,y)}{\partial x}+\mathcal O(h^2)=-\frac{4y}{(4+y^2)^2}+\mathcal O(h^2)$$
it's 2nd order, right? Just as the Laplace equation I approximated in the OP. I mean I've done this a lot of times by now, it's definitely 2nd order right? So why did my programming results show a pretty definitive 1st order approximation?

 

[I like Serena]

Hello again.

I got some help with this since I last posted, but I ran into a problem later on (when calculating all the points on the models when N=M=4,6,8,..,16) which suggested the whole job was done with the 1st level of approximation. So could the programming method (yes, there was programming) be responsible for that? (It was not done with iterations or anything, like some coursemates did it) Or could I have approximated the essential scheme wrong?

Same boundary condition
$$\frac{\partial u(1,y)}{\partial x}=-\frac{4y}{(4+y^2)^2}$$
approximated with central difference.
$$u(1-h,y)=u(1,y)-h\frac{\partial u(1,y)}{\partial x}+\frac{h^2}{2}\frac{\partial^2 u(1,y)}{\partial x^2}+\mathcal O(h^3)$$
$$u(1+h,y)=u(1,y)+h\frac{\partial u(1,y)}{\partial x}+\frac{h^2}{2}\frac{\partial^2 u(1,y)}{\partial x^2}+\mathcal O(h^3)$$
$$\frac{u(1+h,y)-u(1-h,y)}{2h}=\frac{\partial u(1,y)}{\partial x}+\mathcal O(h^2)=-\frac{4y}{(4+y^2)^2}+\mathcal O(h^2)$$
it's 2nd order, right? Just as the Laplace equation I approximated in the OP. I mean I've done this a lot of times by now, it's definitely 2nd order right? So why did my programming results show a pretty definitive 1st order approximation?

Yes. It is second order.
There is 1 assumption though.
That is that $u'''(1 + \theta h, y)$ is bounded for $-1 \le \theta \le 1$.
I mean that its absolute value on the entire interval must be smaller than some (reasonable) constant.
If you hit a singular value, you're in trouble.

Otherwise, I do not have enough information to say anything about where the trouble might be.

 

[I like Serena]

I do see another possibility.

You have introduced fictive nodes for the boundary condition.
Do you solve those fictive nodes for the Laplace equation?
That would be necessary, because otherwise those fictive nodes are left "hanging" and they will generate bad results for the boundary condition.
I think that would match the symptoms.

 

[Thorra]

Well, the fictive nodes are only necessary at the bondaries (upper ones), so yeah, when the laplace equation steered into a boundary point, it transformed the $u_{i+1,j}$ and $u_{i,j+1}$ variables to specific $u_{i,M+1}$ and $ u_{N+1,j}$ fictive points, which have their own special formulas $$u_{i,M+1}=u_{i,M-1}+\frac{2g}{(1+x_i)^2+1}-\frac{4g}{(1+x_i)^2+1}u_{i,M}$$
$$u_{N+1,j}=u_{N-1,j}-\frac{8hy_j}{(4+y^2_j)^2}$$
If you notice I did assume in the end that it's $(1+x_i)+1$ rather than $(1+x_i^2+1)$ cause then the results were off, plus it looks dumb.

 

[I like Serena]

Well, the fictive nodes are only necessary at the bondaries (upper ones), so yeah, when the laplace equation steered into a boundary point, it transformed the $u_{i+1,j}$ and $u_{i,j+1}$ variables to specific $u_{i,M+1}$ and $ u_{N+1,j}$ fictive points, which have their own special formulas $$u_{i,M+1}=u_{i,M-1}+\frac{2g}{(1+x_i)^2+1}-\frac{4g}{(1+x_i)^2+1}u_{i,M}$$
$$u_{N+1,j}=u_{N-1,j}-\frac{8hy_j}{(4+y^2_j)^2}$$
If you notice I did assume in the end that it's $(1+x_i)+1$ rather than $(1+x_i^2+1)$ cause then the results were off, plus it looks dumb.

Huh?
It should be $(1+x_i^2)+1 = 1+x_i^2+1$...
It's derived from $(1+x_i^2) + y^2$, where $y=1$.

 

[Thorra]

But when I do the calculations then it's only close to the analytic results if I change it the way I did.. (Though there's still a strange offset when I go along x-axis with y being constant)

 

[I like Serena]

But when I do the calculations then it's only close to the analytic results if I change it the way I did.. (Though there's still a strange offset when I go along x-axis with y being constant)

Then something else is wrong.
I verified whether the analytic solution matches the Laplace equation and all boundary equations.
It matched (assuming I didn't make a mistake).
The strange offset might be explained if the fictive nodes are not solved using both the Laplace equation and the boundary condition.