Heat equation¶
Shared by Antoine Rideau
On this page you will find how to code using Castor an example as simple as heat diffusion in 1 dimension without heat source with Dirichlet boundary described by the following equation :
The simulation is focuses on the interval \(\left [ x_{min}, x_{max}\right ]=\left [ 0,2 \right ]\) for the space domain which is divided between Nx = 1000
points .
//Parameters
int Nx = 1000;
double xmin = 0;
double xmax = 1;
double tf = 0.1;
double ti = 0;
The space is discretized with dx
steps and the time with dt
steps.
The X
vector stores the space grid
and the vector U0
contains the initial heat distribution
//Discretization
double dx = (xmax - xmin) / (Nx - 1);
double dt = 0.5 * (dx * dx) / d;
matrix<> X = linspace(xmin, xmax, Nx);
matrix<> U0 = pow(sin(X * M_PI), 2);
The second space derivative is approximated by :
Explicit Euler¶
The specifity of the explicit Euler scheme is that the time derivative is calculated with the forward differential quotient approximation
So as to maintain the stability of the simulation, the time step dt
has to respect the following inequality
Finally, the following numerical scheme is obtained :
When coming to the code, \(u_{i}^{n+1}\) is expressed in function of the other dependencies of \(u\) :
double alpha = (dt * d / (dx * dx));
for (double t = 0; t <= tend; t += dt)
{
for (int i = 1; i < Nx - 1; i++)
{
U(i) += alpha * (U(i - 1) - 2 * U(i) + U(i + 1));
}
}
Here you have all the code at once :
#include "castor/matrix.hpp"
#include "castor/graphics.hpp"
using namespace castor;
int main(int argc, char *argv[])
{
//Thermal diffusivity
double d = 1.;
//Parameters
int Nx = 1000;
double xmin = 0;
double xmax = 2;
double tend = 0.1;
double ti = 0;
//Discretization
double dx = (xmax - xmin) / (Nx - 1);
double dt = 0.5 * (dx * dx) / d;
double alpha = (dt * d / (dx * dx));
matrix<> X = linspace(xmin, xmax, Nx);
matrix<> U0 = pow(sin(X * M_PI), 2);
std::cout << "--- Start explicit Euler scheme ---" << endl;
tic();
auto U = U0;
for (double t = 0; t <= tend; t += dt)
{
for (int i = 1; i < Nx - 1; i++)
{
U(i) += alpha * (U(i - 1) - 2 * U(i) + U(i + 1));
}
}
toc();
//Plot
figure fig;
plot(fig, X, U0, {"b-o"}, {"initial"});
plot(fig, X, S, {"g-s"}, {"solution"});
plot(fig, X, U, {"m-x"}, {"explicit"});
drawnow(fig);
return 0;
}
With this code you should get these outputs :
--- Start explicit Euler scheme ---
Elapsed time is 0.213486 seconds.
Implicit Euler¶
The specifity of the implicit Euler scheme is that the time derivative is calculated with the backward differential quotient approximation
dt
.where A is the \(N_{x}\) x \(N_{x}\) tridiagonal matrix
This equation leads to the following linear equation
double alpha = (dt * d / (dx * dx));
matrix<> e = ones(Nx, 1);
smatrix<> A = spdiags(cat(2, cat(2, e, -2 * e), e), colon(-1, 1), Nx, Nx);
auto B = speye(Nx) - alpha * A;
for (double t = 0; t <= tend; t += dt)
{
U = linsolve(B, U);
}
Here you have all the code at once :
#include "castor/matrix.hpp"
#include "castor/graphics.hpp"
#include "castor/linalg.hpp"
using namespace castor;
int main(int argc, char *argv[])
{
//Thermal diffusivity
double d = 1.;
//Parameters
int Nx = 1000;
double xmin = 0;
double xmax = 2;
double tend = 0.1;
double ti = 0;
//Discretization
double dx = (xmax - xmin) / (Nx - 1);
double dt = 5 * (dx * dx) / d;
double alpha = (dt * d / (dx * dx));
matrix<> X = linspace(xmin, xmax, Nx);
matrix<> U0 = pow(sin(X * M_PI), 2);
std::cout << "--- Start implicit Euler scheme ---" << endl;
auto U = transpose(U0);
tic();
matrix<> e = ones(Nx, 1);
smatrix<> A = spdiags(cat(2, cat(2, e, -2 * e), e), colon(-1, 1), Nx, Nx);
auto B = speye(Nx) - alpha * A;
for (double t = 0; t <= tend; t += dt)
{
U = linsolve(B, U);
}
toc();
U = transpose(U);
//Plot
figure fig;
plot(fig, X, U0, {"b-o"}, {"initial"});
plot(fig, X, S, {"g-s"}, {"solution"});
plot(fig, X, U, {"r-+"}, {"implicit"});
drawnow(fig);
return 0;
}
With this code you should get these results :
--- Start implicit Euler scheme ---
Elapsed time is 4.11192 seconds.
References¶
https://www.ljll.math.upmc.fr/ledret/M1English/M1ApproxPDE_Chapter6-2.pdf