Question on heat equation 1D Forward in Time Centered in Space
4 次查看(过去 30 天)
显示 更早的评论
Hi all I have the following question I am trying to solve the PDE with forward time centered in space with the following parameters:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/191039/image.jpeg)
My code so far
function E=Expheat(h,k)
tmax = 1;
t = linspace(0,tmax);
N = 1/h;
M = 1/k;
x = linspace(-pi/2, pi/2);
g(x) = sin*(pi/2*x) + 1/2*sin*(2*pi*x);
alpha = 1;
mu = alpha^2*k/h^2;
%loop for IC
for j = 1 : M
U( 1, j) = 0
end
%loop for
for n = 1: tmax
U(end, j)= exp(pi^2*t/4)
end
for i = 1 : M+1
U(g(x), 0) = (1 - 2*U(n,j))-U(n+1)+U(n-1,j+1)
for j = 2 : N
U(g(x), 0) = (1 - 2*U(n,j))-U(n+1)+U(n-1,j+1)
end
end
I am getting errors and am tying to figure out how to solve this any help would be greatly appreciated.
Thank you
4 个评论
Abraham Boayue
2018-6-6
Here is a code that you may find useful to help solve your problem. I can't really say much about the solution since you did not post the original problem; my code is based on the numerical equations you posted. I would have like to see how the boundary conditions were given instead of the two boundary functions you provided. For example, you may have u(x=0,t) = O1(t), u(x=1,t)= O2(t) or u(0,t) = O1(t), ux(1,t) = O2(t) or something else. If this isn't right, post the original problem as it was given.
clear all
close all
n = 500;
Lx = 1;
dx = Lx/(n-1);
x = 0:dx:Lx;
% 2. Parameters for the t vector
m = 100;
tf = 5;
dt = tf/(m-1);
t = 0:dt:tf;
% 3. Other parameters needed for the solution
% The value of alpha
Fo = 1/4; % a mutiplicative constant that should be < = 1/2
% to insure stability
% Initial and boundary conditions
f = @(x)sin(pi/2*x)+(1/2)*sin(2*pi*x);
g1 = @(t)0;
g2 = @(t)exp(-pi^2/4*t);
u = zeros(n,m);
u(2:n,1) = f(x(2:n)); % Put in the initial condition starting from
% 2 to n-1 since f(0) = 0 and f(N) = 1 for
% N = n-1
u(1,1:m) = g1(t); % The boundary conditions, g1 and g2 at
u(n,1:m) = g2(t); % x = 0 and x = 1
% Implementation of the explicit method
for k = 2:m-1 % Time Loop
for i= 2:n-1 % Space Loop
u(i,k+1) = Fo*(u(i-1,k)+u(i+1,k))+(1-2*Fo)*u(i,k);
end
end
plot(x,u,'--','linewidth',2);
a = ylabel('Pressure');
set(a,'Fontsize',14);
a = xlabel('x');
set(a,'Fontsize',14);
a=title(['Using The Explicit Method - Fo =' num2str(Fo)]);
set(a,'Fontsize',16);
grid;
figure
[X, T] = meshgrid(x,t);
s2 = mesh(X',T',u);
title(['3-D plot of the 1D Heat Equation using the Explicit Method - Fo =' num2str(Fo)])
set(s2,'FaceColor',[0 0 1],'edgecolor',[0 0 0],'LineStyle','--');
a = title('Exact solution of the 1D Diffusivity Equation');
set(a,'fontsize',14);
a = xlabel('x');
set(a,'fontsize',20);
a = ylabel('y');
set(a,'fontsize',20);
a = zlabel('z');
set(a,'fontsize',20);
% disp(u');
采纳的回答
Abraham Boayue
2018-6-7
编辑:Abraham Boayue
2018-6-7
Here is an updated version of the code. I also used matlab pdepe function to validate the results which seem to agree with one another. However, the result obtained from matlab pdepe is more superior than the finite difference method.
clear variables
close all
n = 20;
Lx = 1;
dx = Lx/(n-1);
x = 0:dx:Lx;
% 2. Parameters for the t vector
m = 100;
tf = 1;
dt = tf/(m-1);
t = 0:dt:tf;
% 3. Other parameters needed for the solution
% The value of alpha
Fo = .5; % a mutiplicative constant that should be < = 1/2
% to insure stability
% Initial and boundary conditions
f = @(x)sin(pi/2*x)+(1/2)*sin(2*pi*x);
g1 = @(t)0;
g2 = @(t)exp(-pi^2/4*t);
u = zeros(n,m);
u(2:n,1) = f(x(2:n)); % Put in the initial condition starting from
% 2 to n-1 since f(0) = 0 and f(N) = 1
u(1,:) = g1(t); % The boundary conditions, g1 and g2 at
u(n,:) = g2(t); % x = 0 and x = 1
% Implementation of the explicit method
for k = 2:m-1 % Time Loop
for i= 2:n-1 % Space Loop
u(i,k+1) = Fo*(u(i-1,k)+u(i+1,k))+(1-2*Fo)*u(i,k);
end
end
figure
hold all
for i=1:4:numel(t)
plot(x,u(:,i),'linewidth',1.5,'DisplayName',sprintf('t = %1.2f',t(i)))
end
a = ylabel('Heat');
set(a,'Fontsize',14);
a = xlabel('x');
set(a,'Fontsize',14);
a=title(['Using The Explicit Method - Fo =' num2str(Fo)]);
set(a,'Fontsize',16);
legend('-DynamicLegend','location','bestoutside');
grid;
figure
[X, T] = meshgrid(x,t);
s2 = mesh(X',T',u);
title(['3-D plot of the 1D Heat Equation using the Explicit Method - Fo =' num2str(Fo)])
set(s2,'FaceColor',[0 0 1],'edgecolor',[0 0 0],'LineStyle','--');
a = title('Exact solution of the 1D Diffusivity Equation');
set(a,'fontsize',14);
a = xlabel('x');
set(a,'fontsize',20);
a = ylabel('y');
set(a,'fontsize',20);
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/189073/image.png)
a = zlabel('z');
set(a,'fontsize',20);
% disp(u');
Here is the result from the pdepe function.
function PDEPE_diff
m = 0;
n = 20;
Lx = 1;
dx = Lx/(n-1);
x = 0:dx:Lx;
N = 100;
tf = 1;
dt = tf/(N-1);
t = 0:dt:tf;
sol = pdepe(m,@pdepfun,@icfun,@bcfun,x,t);
u = sol(:,:,1);
figure
hold all % the use of all makes the colors cycle with each plot
for i=1:4:numel(t)
plot(x,sol(i,:),'linewidth',1.5,'DisplayName',sprintf('t = %1.2f',t(i)))
end
title('Numerical solution computed using PDEPE function.')
xlabel('Distance x')
ylabel('Time t')
legend('-DynamicLegend','location','bestoutside');
grid
figure
s2 = surf(x,t,u);
set(s2,'FaceColor',[0 0 1],'edgecolor',[0 0 0],'LineStyle','--');
title('u1(x,t)')
xlabel('Distance x')
ylabel('Time t')
function [g,f,s] = pdepfun(x,t,u,DuDx)
g = 1;
f = DuDx;
s = 0;
function u0 = icfun(x)
u0 = sin(pi/2*x)+(1/2)*sin(2*pi*x);
function [pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t)
pl = ul;
ql = 0;
pr = ur-exp(-pi^2/4*t);
qr = 0;
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/189075/image.png)
更多回答(1 个)
Rauhussaba Rauhy
2022-3-24
How to Build an algorithm for piecewise linear Rayleigh ritz method? _d/dx(p(x)d/dx)+q(x)y=f(x), for x E [0,1] with y(0)=y(1)=0 with piece wise linear function
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!