MATLAB: Hi, I’m trying to solve the heat eq using the explicit and implicit methods and I’m having trouble setting up the initial and boundary conditions.The equation is : du/dt=d^2u/dx^2, initial condition u(x,0)=x, boundary conditions u(0,t)=1 du/dx(1,t)

heat equationnumerical methods

The heat equation I am trying to solve looks like this
du/dt=d^2u/dx^2
initial condition u(x,0)=x
boundary conditions u(0,t)=1 du/dx(1,t)=1
I have attached my work with file name expHeat.m to this question. Thank you

Best Answer

The matlab code posted here is based on the formulation of the explicit method of the finite difference method. I tried to compare the solution to that obtained from using matlab's pdepe solver to ensure that the coding was done correctly. The both results seem to agree, but the solution based on matlab's pdepe solver came out superior in that the mesh points are much finer than the explicit method; which is true since the pdepe function is highly accurate. You can experiment with the parameters to see if you can get a finer mesh points.
The output of the code.
Output of matlab's pdepe solver.
    %%Numerical solution of the Heat equation using the explicit Finite
    % difference method
    clear variables
    close all
    %1 Function parameters
    N  = 500;
    Lx = 50;
    dx = Lx/(N-1);
    x  = 0:dx:Lx;
    alpha = .25;   % To insure stability alpha = dt/dx^2 < = 1/2
    % 2. Time vector
    M = 300;
    tf = 100;     
    dt = tf/(M-1);
    t = 0:dt:tf;
    % 4. Initial and boundary conditions
    f  = @(x) x;   % initial cond. f(x)
    g1 = @(t) 1;     % boundary conditions g1(t) and g2(t)
    g2 = @(t) 1;
    % 5. Inilization of the heat equation
    u      =  zeros(N,M);
    u(:,1) =  f(x);
    u(1,:) =  g1(t);
    % 6. Implementation of the explicit method
    for j= 1:M-1       % Time Loop 
         for i= 2:N-1  % Space Loop
           u(i,j+1) = alpha*(u(i-1,j))+(1-2*alpha)*u(i,j) + alpha*u(i+1,j);
         end
         % Vectorize the inner for loop with this line to gain some speed.
         % u(2:N-1,j+1) = alpha*(u(1:N-2,j))+(1-2*alpha)*u(2:N-1,j) + alpha*u(3:N,j);
     % Insert boundary conditions for i = 1 and i = N here.
     u(2,j+1) =  alpha*u(1,j) + (1-2*alpha)*u(2,j) + alpha*u(3,j);
     u(N,j+1) =  2*alpha*u(N-1,j)+(1-2*alpha)*u(N,j)+ 2*alpha*dx*g2(t);
    end
    % Plot results
    figure
    plot(x,u(1:end,1:30),'linewidth',2);
    a = ylabel('Temperature');
    set(a,'Fontsize',14);
    a = xlabel('x');
    set(a,'Fontsize',14);
    a=title(['Using The Explicit Method - alpha =' num2str(alpha)]);
    legend('Explicit soln')
    set(a,'Fontsize',16);
    xlim([0 1]);
    grid;
    figure
    [X, T] = meshgrid(x,t);
    s2 = mesh(X',T',u);
    title(['3-D plot of the 1D Heat Eq. using the Explicit Method - alpha = ' num2str(alpha)])
    set(s2,'FaceColor',[1 0 1],'edgecolor',[0 0 0],'LineStyle','--');
    a = title('Explicit finite difference 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);
    xlim([0 2])
    zlim([0 2.5])
Related Question