Numerically solving a diffusion equation with a piecewise initial condition

22 次查看（过去 30 天）
Bas123
Bas123 2023-4-25
评论： Alexander 2023-12-2
Ran in:
Consider the diffusion equation given by with initial conditions and boundary conditions .
I wish to numerically solve this equation by the finite difference formula where with ∆x = 0.1 and ∆t = 0.01. The exact solution is given by
Here is the code that I have worked out so far.
clear all
D = 1/4;
dx = 0.1;
dt = 0.01;
x = 0:dx:1;
t = 0:dt:1;
Nx = length(x);
Nt = length(t);
u = zeros(Nx, Nt);
u(x<=1/2,1) = x(x<=1/2);
u(x>1/2,1) = 1 - x(x>1/2);
u(1,:) = 0;
u(Nx,:) = 0;
for j = 1:Nt-1
for i = 2:Nx-1
u(i,j+1) = D*(dt/dx^2)*(u(i+1,j) + u(i-1,j)) + (1-2*D*(dt/dx^2))*(u(i,j));
end
end
% Compute the analytical solution V(x,t)
k = 1:100;
V = zeros(length(x), length(t));
for j = 1:Nt
for i = 1:Nx
V(i,j) = sum((4./(k.*pi).^2).*sin(k.*pi/2).* sin(k.*pi.*x(i)).*exp(-D.*t(j).*(k.*pi).^2));
end
end
% Plot the numerical and exact solutions
figure;
[X,T] = meshgrid(x,t);
surf(X,T,u');
xlabel('x');
ylabel('t');
zlabel('u(x,t)');
title('Numerical solution');
figure;
surf(X,T,V');
xlabel('x');
ylabel('t');
zlabel('V(x,t)');
title('Exact solution');
However, there seems to be something wrong with the code. Any support would be greatly appreciated.
  3 个评论
显示 1更早的评论
Bas123
Bas123 2023-4-25
Hi Bill, thanks for your message. I am supposed to use an explicit finite difference method for the problem.
Alexander
Alexander 2023-12-2
Hi! This was really helpful code for a HW assignment I'm working on. My graphs are a little wacky as well, I am using different differencing schemes, but I think the problem may lie in your parameters. I think by finding the correct dx and dt you may be able to generate smoother graphs.
However I am also pretty new to this and not sure exactly how to find the best spatial and time step sizes. Good luck!

请先登录，再进行评论。

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

采纳的回答

Torsten
Torsten 2023-4-25
编辑：Torsten 2023-4-25
I changed
V(i,j) = sum((4./(k.*pi).^2).*sin(k.*pi/2).* sin(k.*pi.*x(j)).*exp(-0.5.*t(i).*(k.*pi).^2));
to
V(i,j) = sum((4./(k.*pi).^2).*sin(k.*pi/2).* sin(k.*pi.*x(i)).*exp(-D.*t(j).*(k.*pi).^2));
in your code and the results look at least similar.
You should check the results in detail, i.e. plots of u over x for some times t or plots of u over t for some positions x. A surface plot looks fine, but is not suited to find errors or unprecise results.

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Eigenvalue Problems 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by