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

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

[Jan]

Please post the error messages.

[Janvier Solaris]

>> E=Expheat(10,5) Error using sin Not enough input arguments.

Error in Expheat (line 7) g(x) = sin*(pi/2*x) + 1/2*sin*(2*pi*x);

[Abraham Boayue]

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]

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

Answer 1

[Abraham Boayue]

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;

[Janvier Solaris]

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]

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]

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