0
$\begingroup$

I am trying to create a function that solves one BTCS step of the following PDE:

$\frac{\partial q}{\partial t} = k\frac{\partial^{2}q}{\partial x^{2}}$

over a domain of: $0\leq x \leq \pi$

with initial condition given by the following square wave:

$q(x,0) = \left\{ \begin{array}{ll} 1, & \frac{2\pi}{5}<x<\frac{3\pi}{5}\\ 0, & otherwise \end{array} \right.$

With boundary conditions:

$q(0,t) = 0$

$\frac{\partial q}{\partial x}(\pi,t) = 0$

and $k = 0.05$

Here is my code:

function [ qnew ] = btcs_step( q,dt,dx,k )
    %INPUTS:
    % q: column vector of solution values from previous timestep
    % dt: timestep
    % dx: x-step
    % k: diffusion equation coefficient 

    % generate A matrix
    A_size = length(q);
    A = gallery('tridiag', A_size, -k/dx^2, 1/dt + (2*k)/dx^2, -k/dx^2);

    % generate b matrix
    %b = ones(length(q),1).*(1/dt).*q';
    b = q'/dt;

    %impose boundary conditions
        % Dirichlet conidtions
    b(1) = 0;
    b(end) = b(end-1);

        % Modifications to coeff matrix
    A(1,1) = 1;
    A(1,2) = 0;
    A(end,end-1) = 0;
    A(end,end) = 1;

    % solve matrix
    qnew = A\b;
end

The second Dirichlet condition is a finite difference approximation to the derivative

I am calling this function and plotting results for various timesteps at various times with the following code:

% parameters
k = 0.05;
N = 100;

dx = pi/N;
timesteps = [0.005, 0.5*dx^2/k, 0.1];
x_vals = 0:dx:pi;

%initial condiiton
q_init = (2*pi/5)<x_vals & x_vals<(3*pi/5);

%stuff for plotting
plot_times = [0,2.5,5,10,20,40,80];
legend_text = strcat({'t = '},string(plot_times));


for dt = timesteps
    % create data structure
    tvals = 0:dt:80; nt = length(tvals);
    q = zeros(nt, N+1);
    q(1,:) = q_init;
    for t = 2:nt
        q(t,:) = btcs_step( q(t-1,:), dt, dx, k);
    end
    figure;
    hold on;
    for pt = plot_times
        pti = floor(pt/dt)+1;
        plot(x_vals,q(pti,:));
    end
    title(sprintf('FTCS Solution of 1-D diffusion equation\nk = 0.05\nTimestep: %g',dt));
    legend(legend_text);
    xticks([0,pi/4,pi/2,3*pi/4,pi]);
    xticklabels({'0','\pi/4','\pi/2','3\pi/4','\pi'});
    axis([0,pi,-inf,inf])
    hold off
end

Here are the three plots:

plots

There should be 7 curves on each plot...

$\endgroup$
  • 1
    $\begingroup$ There are possibilities of errors at a mathematical level (using the wrong equations) and at the programming level (incorrect implementation of the equations.) You should add a discussion of the equations that you're using to your question. $\endgroup$ – Brian Borchers May 4 '18 at 15:32

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Browse other questions tagged or ask your own question.