Numerical integration of PDEs: How do you satisfy boundary conditions

[davidbenari]

Suppose we are solving a diffusion equation.

$\frac{\partial}{\partial t} T = k\frac{\partial^2}{\partial x^2} T$

On the domain $0 < x < L$

Subject to the conditions

$T(x,0) = f(x)$ and $T = 0$ at the end points.

My question is:

Suppose we solve this with some integration scheme (forward time centered space), such that we use a formula like

$T_{i}^{n+1} = \frac{k\tau}{h^2} (T_{i+1}^{n} - 2 T_{i}^{n} + T_{i-1}^n)+T_{i}^n$

how is it that we are ensuring that $T=0$ at the end points in the future?

Is a partial differential equation completely determined by its initial configuration (including its initial boundary)?

How are we respecting the boundary conditions when we do such integration of the pde ?

 

[Twigg]

That finite difference scheme you supplied is a central difference scheme, meaning it applies exactly as written only for spatial mesh points (x points) between the end points (exclusively). Its easiest to supply the boundary conditions if you treat that recurrence relation as a matrix (or an excel formula). For consistency, let's say each column holds data varying over time at constant position and each row contains data sampled over space (x coordinate) at a fixed time. Then your initial condition ( t = 0 ) fills all the data in the top row ( n = 0 ). This isn't enough information to fully determine the unknown data. If you put that recursion formula into a rectangular table in excel and filled the top row only with givens, it wouldn't find any data when it looks left ( i - 1 ) and right ( i + 1 ). It would only see the data above it ( n - 1 ).
The end points correspond to the leftmost ( i = i_min, x = 0 ) and rightmost ( i = i_max, x = L) columns of the array. What I've done in the past and worked in MATLAB is to overwrite these side columns with the fixed boundary value data (in your case, T = 0 for all entries in both the first and last columns). Then you have 1 temporal constraint and 2 spatial constraints for a difference scheme that is 1st order in time and 2nd order in space. It's fully defined.

 

[davidbenari]

Hmm. So you only overwrite the boundary points?

 

[Lord Crc]

You're not supposed to ensure the boundary conditions hold as such, the boundary conditions tell you what the value should be.

Think of a guitar string. No amount of plucking will move the endpoints of the string. They're fixed in place.

 

[davidbenari]

You're not supposed to ensure the boundary conditions hold as such,

When you solve these problems analytically you basically force your solution to be zero at the end points. Think of when, when applying separation of variables and you get some sine function as a solution, you construct the sine's argument to ensure that its zero on the endpoints.

When we solve PDEs numerically, I don't see how were forcing our boundaries to hold a specific condition. It seems logical to me to just overwrite the boundaries in each iteration step though.

 

[Lord Crc]

Disclaimer: it has been a few years since I left uni, so maybe a bit rusty...

When you solve these problems analytically you basically force your solution to be zero at the end points. Think of when, when applying separation of variables and you get some sine function as a solution, you construct the sine's argument to ensure that its zero on the endpoints.

Yes that's what I mean. You assume that nothing your solution does affects that value, and so you use the value as-is at the boundary. Thus your solution does not ensure the boundary conditions hold, it assumes they do.

When we solve PDEs numerically, I don't see how were forcing our boundaries to hold a specific condition.

Let's write out explicitly how to compute T^{n+1}_1 given your above formula:T^{n+1}_1 = \frac{k\tau}{h^2}\left(T^n_2 - 2T^n_1 + T^n_0\right) + T^n_1

If n &gt; 1, what value will you use for T^n_0? You can't compute it using your formula, because then you'd need to know the value of T^n_{-1} which is outside the domain and thus not defined.

But, you have a boundary condition saying that T^n_0 = 0. So let's use that: T^{n+1}_1 = \frac{k\tau}{h^2}\left(T^n_2 - 2T^n_1 + 0\right) + T^n_1T^{n+1}_1 = \frac{k\tau}{h^2}\left(T^n_2 - 2T^n_1\right) + T^n_1
Now you've forced the value at T^n_0 to be zero, and used that to solve the PDE, instead of ensuring it winds up as zero afterwards.

Of course, in this case you could have just overwritten the boundary values in your array holding T in your program, as Twigg suggested, and used your regular formula for all internal points. However in some cases you have more complicated boundary conditions where it makes more sense to change the expression used to evaluate the points near the boundaries, like I did above for T^{n+1}_1.