Solving the heat equation using finite difference with tridiagonal matrix

I have the following code to solve the heat equation via a tridiagonal matrix but have no clue how to continue.

clear
nx=10;
a=0; b=1;
Du=1; dt = 0.1;
alpha=1; beta=1; %BCs
p=0; q=0;
x=linspace(a,b,nx+1); dx=(b-a)/nx; %forms the grid and computes the step-size
u(1)=alpha; u(nx+1)=beta; %sets the boundary conditions
A=tridiag(1/(dx*dx)-p/(2*dx),q-2/(dx*dx),1/(dx*dx)+p/(2*dx),nx-1);
% A is the tri-diagonal matrix containing the LHS of the finite diff equations
r=functionr(x)';
% r gives the RHS function.
% The Dirichlet bc must be incorporated in the 1st and
% last equation of the system as follows:
r(2)=r(2)-(1/(dx*dx)-p/(2*dx))*alpha;
r(nx)=r(nx)-(1/(dx*dx)+p/(2*dx))*beta;
% u(1) and u(n+1) given by BC,
% values in between computed from solving the linear system as follows:
% Euler's Method
h = (Du*dt)/dx^2;
for j=1:nx-1 % loop for Euler's method
    u(j+1) = u(j) + h*A*u(j);
end
plot(x,u)

the tridiag function is as follows

function T = tridiag(a, b, c, n)
% tridiag Tridiagonal matrix.
% T = tridiag(a, b, c, n) returns an n by n matrix that has
% a, b, c as the subdiagonal, main diagonal, and superdiagonal
% entries in the matrix.
T = b*diag(ones(n,1))+c*diag(ones(n-1,1),1)+a*diag(ones(n-1,1),-1);

the functionr is as follows

function r = functionr(x)
r = x*0;

all help much appreciated

[Torsten]

We do not know your problem, not to mention how to continue.

Answer 1

[Rohit Kulkarni]

Hi Joshua,

In my understaing it looks like you're trying to solve a differential equation using a finite difference method, but there seems to be some issues in your code:

Issues and Improvements:

1. Initialization of `u`:You need to initialize `u` for all grid points, not just the boundaries. Initially, you can set it to zeros or any initial condition if given.

2. Matrix-vector multiplication: In the Euler method loop, the way you're trying to update `u(j+1)` is incorrect. You should perform a matrix-vector multiplication for the whole vector `u` (excluding the boundary values) at once, not element by element.

3. Boundary conditions handling: The boundary conditions are correctly set at the beginning, but when updating `u`, you should exclude the boundaries since they are constant (for Dirichlet conditions).

4. Function functionr:The function `functionr` returns a zero vector, which means it doesn't contribute to the solution. If there's supposed to be a source term or any non-zero right-hand side, it should be defined here.

Improved Code:

clear
nx = 10;
a = 0; b = 1;
Du = 1; dt = 0.1;
alpha = 1; beta = 1; % BCs
p = 0; q = 0;
x = linspace(a, b, nx+1); dx = (b-a)/nx; % forms the grid and computes the step-size
u = zeros(nx+1,1); % Initialize u for all grid points
u(1) = alpha; u(nx+1) = beta; % sets the boundary conditions

% Define A matrix
A = tridiag(1/(dx*dx)-p/(2*dx), q-2/(dx*dx), 1/(dx*dx)+p/(2*dx), nx-1);

% Define r vector
r = functionr(x(2:end-1))'; % Exclude the boundary points in x

% Adjust r for boundary conditions
r(1) = r(1) - (1/(dx*dx)-p/(2*dx))*alpha;
r(end) = r(end) - (1/(dx*dx)+p/(2*dx))*beta;

% Time-stepping with Euler's method
h = (Du*dt)/dx^2;
for j = 1:nx-1
    u(2:nx) = u(2:nx) + h*(A*u(2:nx) + r); % Update interior points only
end

plot(x, u)
xlabel('x')
ylabel('u')
title('Solution of the Heat Equation')

Hope this resolves your query!

Thanks