1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
|
% Numerical Analysis
% Group 10: Forbelets
% Imperial College, March 2017
## Grahps Benedikt
{width=80%}
{width=80%}
{width=80%}
{width=80%}
{width=80%}
{width=80%}
{width=80%}
{width=80%}
{width=80%}
{width=80%}
{width=80%}
{width=80%}
{width=80%}
{width=80%}
{width=80%}
# Finite Differences for PDE
## Comments on implementation:
We first create a matrix to implement our finite difference method on. The two dimensions of the matrix and time and displacement.
The first row of the matrix (in other words at time t = 0) is populated with the initial conditions.
The rest of the matrix is populated in a loop from the initial conditions and the boundary conditions.
In the loop, for each value of m, the time iterate of the matrix, we set the first and last possible displacement index's to our desired boundary conditions. We then we iterate trough each displacement step and apply the central algorithm to it. The full implementation of this loop can be seen in the included code in ex4.m.
$$ U_{j}^{m+1} = v U_{j-1}^{m} + (1-2v) U^m_j + v U^{m}_{j+1} $$
$$ = U^m_j + v ( U_{j}^{m+1} - 2 U^m_j + U^{m}_{j+1}) $$
Where $$ v=\frac{k}{h^2} $$
Writing this in matlab as a two dimensional matrix is very similar.
```matlab
u(m+1,j) = u(m,j) + ((k/(h^2))*(u(m,j+1) - 2*u(m,j) + u(m,j-1)));
```
Physicaly this means that if that if the sum of the temperature of the two adjacent points is higher is higher than two times the temperature of the given point, the temperature of the selected point will go up after one time step, while if the two adjacent points sum up to less than twice our current temperature, our next temperature will decrease.
As time tends to infinity the temperature will tend towards the boundary conditions as they will set the temperature of their neghbours which will propagate trough the curve.
## Choosing h and k
Choosing h (the displacement step) and the appropirate k (time step) is not a trivial task. Picking a very small time step is computationally expensive (requires a big matrix and a computer with more memory), and we generally want to compare the temperature distribution at larger time increments (as shown in the instruction sheet).
The Von Neumann stability anlaysis shows that in order to produce a stable results it is required:
$$ \frac{k}{h^2} \leq \frac{1}{2} $$
Since we want to produce long simulations with a small distance step, but don't want to see extremely small steps in the graph, we set the time step to satisfy the equality given a value for h. Even then, the steps are hardly small so when graphing we only a fraction of the time steps from the whole range. In our implementation we print 10 equaliy spaced lines from the time range.
### Tent Function
$$
y_0(x) =
\left\{
\begin{array}{ll}
2x & \mbox{for } x \in [0,0.5] \\
2-2x & \mbox{for } x \in [0.5,1]
\end{array}
\right.
$$
{width=80%}
In figure \ref{tent} we can see the tent function applied as the initial condition to a one dimensional rod of length 1. The coloured lines represent the approximated temperature at times increments of 0.1, all the way up to 1.0. Because of the boundary conditions the end points are bound to zero degrees and as a result the temperature of the rod decreases. As **t** approaches infinity the temperate of the rod will approach zero, where the center point which started as the highest temperature and is furthest away from the boundaries maintains the highest temperature of all points.
### Sinusoidal Function
$$ y_0(x) = \sin{2\pi x} $$
{width=80%}
In figure \ref{sina} we perform the finite difference method on to a sinusoid. The result yet again is decreasing temperature as the boundaries are zero. The temperature decreases such that the resulting distribution continues to be sinusoidal, only with a different amplitude. The left and right side fo the rod have identical but opposite amplitudes, and due to this equilibrium the center point at 0.5 displacement stays at 0 temperature.
$$ y_0(x) = |\sin{2\pi x}| $$
{width=80%}
In figure \ref{sinb} we use the same initial condition function, only this time we take the absolute value of it so what was negative before is now positive. It is interesting to note that the resulting distribution at $t \neq 0$ is very different from the absolute value of that in figure \ref{sinb}. Point 0.5 no longer sits at 0 and is instead heated due to neighboring elements until it becomes the hottest point.
### Square Function
$$
y_0(x) =
\left\{
\begin{array}{ll}
0 & \mbox{for } x \in [0,0.25] \cup [0.75,1] \\
1 & \mbox{for } x \in [0.25,0.75]
\end{array}
\right.
$$
{width=80%}
Applying the "square" function in figure \ref{square} produces interesting results. Where as in figure \ref{sinb} the center heated up while the corners cooled down the opposite happens initially, as the curve straightens out. This happens because of the extreme temperature difference that we introduced.
### Double Tent function
{width=80%}
Making two symmetric tent functions with oposite signs like in figure \ref{dtent},shows a result like that in figure \ref{tent}, with an artificial boundary of zero in the centre due to the symmetry.
## Varying boundary conditions
{width=80%}
It is possible to give the algorithm an initial condition which does not match the boundary conditions. In figure \ref{extra1} we can see the initial conditions in blue which attempt to set the corners of the rod to -1 and 1 respectively. As the boundaries are forced to to zero they quickly pull the graph away from the original initial conditions and towards themselves. As **t** approaches infinity \ref{extra1} will become a straight line between the two boundaries. Again we have an equilibrium and a constant temperature of zero in the middle.
### Non Zero boundary conditions
{width=80%}
In \ref{extra2} we use the same initial conditions as in figure \ref{sina}, however this time we set the boundary conditions to be 0.5 and 1. Yet again this forces the temperature curve towards them. The equilibrium at 0.5 displacement is no longer maintained as as the boundaries are different and non-zero. As **t** approaches infinity we will get a straight line between 0.5 and 1 temperature.
### Time varying boundary conditions
{width=80%}
Finally in every loop itteration we can change the boundary condition based on the elapsed time. We can see the effect of doing so on the original tent function in figure \ref{extra3}. The boundaries move from 0 degrees to 0.5 degrees at a rate of 5 degrees per second. The result is the function bending back on itself as the boundaries start moving up.
|
