solve heat conduction equation pde in matlab using implicit euler for time and central differencing for space in matlab

The Heat conduction equation in one dimension is given as:

d/dt(u(x,t)) = alpha(d^2/dx^2)(u(x,t))

where u(x,t) is the temperature at position x and time t, and alpha is the thermal diffusivity of the material.

To solve this equation using central differencing for space and implicit Euler for time in Matlab, we need to discretize the space and time domains.

Let nx be the number of space points and nt be the number of time points.

We can discretize the space domain with a spacing of delta_x as:

x_i = (i-1)*delta_x, where i = 1,2,3,...,nx

Similarly, we can discretize the time domain with a time step delta_t as:

t_j = (j-1)*delta_t, where j = 1,2,3,...,nt

Let u(i,j) be the numerical solution at position i and time j.

Using central differencing for space, we can approximate the second derivative as:

(d^2/dx^2)(u(i,j)) = (u(i-1,j) - 2*u(i,j) + u(i+1,j))/(delta_x^2)

Using implicit Euler for time, we can approximate the time derivative as:

(d/dt)(u(i,j)) = (u(i,j) - u(i,j-1))/delta_t

Substituting these approximations in the heat equation, we get:

(u(i,j) - u(i,j-1))/delta_t = alpha*(u(i-1,j) - 2*u(i,j) + u(i+1,j))/(delta_x^2)

Solving for u(i,j) at time step j, we get:

u(i,j) = u(i,j-1) + (alpha*delta_t/delta_x^2)*(u(i-1,j) - 2*u(i,j) + u(i+1,j))

This is an implicit equation in u(i,j) as u(i,j) appears on both sides of the equation. Hence, we need to solve a system of equations to get the numerical solution at time step j.

Let A be the coefficient matrix and b be the right-hand side vector of the system of equations.

We can write the equation for u(i,j) in a matrix form as:

A(i,i-1) = -alpha*delta_t/(delta_x^2), A(i,i) = 1 + 2*alpha*delta_t/(delta_x^2), A(i,i+1) = -alpha*delta_t/(delta_x^2), b(i) = u(i,j-1)

We need to solve for u(:,j) using the matrix equation:

A*u(:,j) = b

The initial condition and boundary conditions need to be specified to get a unique numerical solution.

Let u0(x) be the initial temperature distribution and ul(t) and ur(t) be the temperatures at the left and right boundaries at time t, respectively.

We can initialize u(i,1) as u0(x_i).

The boundary conditions can be approximated using forward and backward finite differences as:

u(1,j+1) = ul(t_j+1) and u(nx,j+1) = ur(t_j+1)

The Matlab code to solve the heat equation using implicit Euler and central differencing can be written as:

main.m
% Heat conduction equation using implicit Euler and central differencing

% Parameters
alpha = 1; % Thermal diffusivity
L = 1; % Spatial domain
T = 1; % Time domain

% Spatial discretization
nx = 50; % Number of points
delta_x = L/nx; % Spacing
x = linspace(0, L, nx)'; % Grid

% Time discretization
nt = 100; % Number of time steps
delta_t = T/nt; % Time step
t = linspace(0, T, nt+1)'; % Grid

% Boundary conditions
ul = 0; % Left temperature
ur = 0; % Right temperature

% Initial condition
u0 = sin(pi*x/L); % Initial temperature distribution

% Numerical solution
u = zeros(nx, nt+1); % Solution matrix
u(:,1) = u0; % Initial condition

% Coefficient matrix
r = alpha*delta_t/delta_x^2;
A = diag(1+2*r*ones(nx,1)) + diag(-r*ones(nx-1,1), 1) + diag(-r*ones(nx-1,1), -1);

% Time integration
for j = 1:nt
    b = u(:,j);
    b(1) = ul(t(j+1));
    b(nx) = ur(t(j+1));
    u(:,j+1) = A\b;
end

% Plot solution
surf(x, t, u')
xlabel('x')
ylabel('t')
zlabel('u(x,t)')
968 chars
46 lines

This code discretizes the space and time domains and solves for u(i,j) using an implicit scheme. It also plots the solution as a surface plot.

Note that the stability of this scheme depends on the time step and the size of the spatial discretization. To ensure stability, the time step should be small enough and the spatial discretization should be fine enough.