Question on heat equation 1D Forward in Time Centered in Space

Question

Janvier Solaris 2018-6-5

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/404194-question-on-heat-equation-1d-forward-in-time-centered-in-space

回答： Rauhussaba Rauhy 2022-3-24

Hi all I have the following question I am trying to solve the PDE with forward time centered in space with the following parameters:

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 个评论
显示 2更早的评论隐藏 2更早的评论

Abraham Boayue 2018-6-6

在 MATLAB Online 中打开

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');

Janvier Solaris 2018-6-6

Hi Abraham, Thank You for your reply, here is the original Problem

请先登录，再进行评论。

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

Answer 1

Abraham Boayue 2018-6-7

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/404194-question-on-heat-equation-1d-forward-in-time-centered-in-space#answer_323676

编辑：Abraham Boayue 2018-6-7

在 MATLAB Online 中打开

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);

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;

2 个评论
显示无隐藏无

Janvier Solaris 2018-6-11

Thank you Abraham for providing the pdpe so as to compare, I updated my code with your input and was able to clear my errors Thank You

Harsha 2020-11-10

编辑：Harsha 2020-11-11

9A30E1AF-E2B8-4DEE-8567-6AC738DFBFF1.jpeg

Hi,

I have the same exact question but rather I want to get the results as tabulated values as t versus i or x like the attached pic. Could someone help please. With boundry condtions u(0,0.01) and u(1,0.01) and intial condition f(x)=100

Thanks

请先登录，再进行评论。

Answer 2

Rauhussaba Rauhy 2022-3-24

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/404194-question-on-heat-equation-1d-forward-in-time-centered-in-space#answer_925859

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