How do write the code (or suggestions on where I should begin) to plot every 10th iteration?

Hello,

I am using a numerical scheme (crank-Nicholson scheme) to solve a PDE as part of an assignment. To see what I am exactly doing, I copy and pasted the entire code.

So far I have it as a function and it plots every iteration. I don't think it stores it. Could plotting also be making my function slower?

Any suggestions for improvement would be welcome. :)

This is an diffusion-advection-reaction problem.

The equation in question is the following.

dC/dt= D*[d^2*C/dx^2]- ad*[dC/dx]- R(x)

~Cat

************************************************************

function [pp ] = CrankNick3
%Solves the advection-diffusion-reaction PDE using the crank-nicholson
%scheme
% instead of using the command AA\BB, the conjugate graduent method is used
% to speed up computation
%the conjugate gradient code was obtained from MathWorks
xl=0;xr=1;%boundaries in space
yb=0;yt=1; %boundaries in time
M=50; %number of points in space
N=100; %number of points in time
dx=(xr-xl)/M;dt=(yt-yb)/N;            % step sizes
m=M-1; n=N;
x=xl:dx:xr;tt=yb:dt:yt;
xx=x(1:m);

G0=10; %uM
%coefficients 
global D W k b tol
D=8e-8; %diffusion
W=-1*.1e-2;%advection
tol=1e-5;

%Rxn should be a m*1 vector. The rxn coefficient of the PDE varies in space
k=.15; %Rcox from Burdige
b=.1; %the beta reaction term
Rxn=k*exp(-b*xx'); 
%Coefficients used to generate matrix AA and BB
global alpha beta gama
alpha=dt*D/(2*dx);
beta=-dt*W/(4*dx);
gama=-dt*Rxn./2;

global c1 c2 c3 c4 c5 c6
c1=-alpha-beta;
c2=1+2*alpha-gama(1:m);
c3=-alpha+beta;
c4=alpha+beta;
c5=1-2*alpha+gama(1:m); 
% c5=gama(1:m);
c6=alpha-beta;

%Generate AA matrix
% define tridiagonal matrix AA
AA=diag(c2.*ones(m,1))+diag(c3.*ones(m-1,1),1);
AA=AA+diag(c1*ones(m-1,1),-1); 
%the next few lines define the boundary conditions
AA(1,:)=0; %the first row of any colum must be zero
AA(1,1)=1; %the first row, first column value must be one
AA(m,:)=0; %the last row, any colum must be zero
AA(m,m)=1; %the last row, last column must be one

%form initial matrix B0 with some concentration profile
% B0=G0;
Bb=zeros(m,1); %this forms the initial b matix
Bb(1,1)=(G0); %this defines the initial conditions
w0=AA\Bb; 

BB=zeros(m,1);

BB=Bb; %for the first iteration of BB, let it be initial condition vector Bb
w=w0; %for the first iteration of w, let it be initial condition vector w0
tic
for j=1:n,
    for p=2:m,
        BB(p)=(c4*w(p))+(c5(p).*w(p))-(c6*w(p));

            BB;
            gcf; 
            xp=x';

             fignodiff=plot(xp(1:49,1),w);hold on;

            xlabel('depth');ylabel('c');
            drawnow;
        end
    w=conjgrad(AA,BB,tol);

end
toc
end

[Jan]

Please format your code as explained in the "Markup help" link on this page.

Answer 1

[per isakson]

I red the first two first sentences of your question.

Plot makes is slower, probably much slower.

    ii=ii+1;
    if ii == 10
      set( ..., 'Xdata', ... , 'Ydata', ... )
      ii = 0;
    end

Just resetting the data is faster than plot(...).

--- EDIT ---

It is an investment to learn to use Handle Graphic of Matlab.

Firstly, you need understand how Matlab graphic displays are composed by Handle Graphic Objects. Each graphic object has a handle. (Don't confuse with objects and the handle class of the OOP-system.) The graphic objects have properties - a lot of them. In the online help you will find the entry "Handle Graphics Property Browser", which provides an overview.

.......

The details are in the documentation. You need to study some examples and do some experiments.

    set( line_handle, 'XData', new_x_data, 'YData', new_y_data )

will modify a line object in an existing diagram.

It cannot be done by a few tips and trial and error. You need to invest some time with the documentation (, which I think is good - both spending time and the documentation).

"where I should begin". Run: Demos / Graphics / Line Plotting , which you find under Help / Contents. There you may change values of properties and the corresponding code is shown in a mini command window.

[Cathleen]

I went to the help on set...but I am not sure what set does.

Would you mind explaining what set does/how does it work?

Sorry about asking more questions...

Answer 2

[per isakson]

Now I have run your code with

    profile on

I replaced the conjgrad-line by "w=w+rand" and moved tic&toc inside the outer loop.

I would have done something like this (very old habit):

    fig_handle = figure( ... )
    ax_handle  = axes( 'parent', fig_handle, ... )
    line_handle= line( 'parent', ax_handle, ... )

    ii = 0;
    for jj=1:n,
        for p=2:m,
            ...
            ...
            if ii >= 10
                set( line_handle, 'Ydata', w  )
                drawnow
                ii = 0;
            end 
            ii = ii + 1;             
            ...
        end
    end

The property, DoubbleBuffer, of the figure and EraseMode of the line was once important to simple animations. However, with todays computers I don't know.

You want to add a new lines to the diagram, not "replace" the existing line? There will be many (5000) line objects, which are redrawn every time.

[Cathleen]

Thank you. I tried this but things didn't work out. Instead I saved the results in a "growing" matrix (augments the result after each itr) and plotted it afterwards. Thank you for your help!