Solving the heat equation using finite difference with tridiagonal matrix

30 次查看（过去 30 天）
Joshua
Joshua 2024-3-9
回答： Rohit Kulkarni 2024-3-20
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

请先登录，再回答此问题。

回答（1 个）

Rohit Kulkarni
Rohit Kulkarni 2024-3-20
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

类别

在 Help Center 和 File Exchange 中查找有关 Partial Differential Equation Toolbox 的更多信息

标签

产品


版本

R2022b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by