Help with CODE please

Hi,

I am trying to model this diffusion equation but can't seem to get the code to work. Was hoping you guy can give me a hand.

The equation is as follows:

C(r,t)=(Co/2Dt)exp(-r^2/(4Dt))

Where, 
D=6E-10

I want the user to be able to pick the time and the initial concentration (Co). Then plot C vr r. 

This is what I have so far:

%{
SPECIFY parameters
%}

t=input('Enter time in seconds:');            % Request user input for time
Co=input('Enter initial NO Concentration:');  % Request user input for Co

D=6E-10;        % Diffusivity Coefficient, m2s-1

i=1;            % Operator looking into the pre-located y array 

y=(1:1:2000);
%{
RUNNING the simulation
%}

for x=(0:0.0005:0.5)    %For x values between 0 and 50cm

      y(i)=(Co/(2*D*t))*exp(-((x^2)/(4*D*t)));

     i=i+1;
  end
%{
PLOTTING results
%}

plot(y);
xlabel('Height (m)');
ylabel('NO Concentration (M)');

Any help would be very much appreciated. Thnx

[Geoff Hayes]

Mo - please format the above code so that it is readable. Highlight the code portions and press the {}Code button.

As well, please indicate what in particular is causing a problem.

A couple of things I did notice: the code will fail when you try to add the labels to the x and y axes of your plots - the commands are xlabel and ylabel respectively. Consider creating an x vector using linspace to size x to be the same as y (for values from 0 to 50 centimetres). Then iterate over the length of x and plug x(i) into the equation. With an actual vector x, you can use this in your plot command as plot(x,y). Please check the documentation of this command (type help plot) to see why this is different than what you have.

[Image Analyst]

Geoff, put it in an Answer so it can be accepted, and you can get credit.

[Mo]

Geoff Hayes - I know how to create the x vector using the linespace to size x but not sure how I can iterate over the length of x and plug x(i) into the equation. If you could help me with that, I would appreciate it.

This is my so not successful attempt:

%{
SPECIFY parameters
%}

t=input('Enter time in seconds:');          % Request user input for time
Co=input('Enter initial NO Concentration:');% Request user input for Co

D=6E-10;        % Diffusivity Coefficient, m2s-1

i=1;            % Operator looking into the pre-located y array
%{
RUNNING the simulation
%}

for x=linspace(0,0.1,0.0005);    % x values between 0m and 0.1m

      y(i)=(Co/(2*D*t))*exp(-((x^2)/(4*D*t)));

     i=i+1;
  end
%{
PLOTTING results
%}

plot(x,y);       
xlabel('Height (m)');
ylabel('NO Concentration (M)');

[Geoff Hayes]

Mo - originally you had the y vector set to 1x2000 element array (no need to initialize it from 1 to 2000 since the values are going to be overwritten). So you can try something like:

n = 2000;       % if 2000 is too large, then you can just change n
y = zeros(1,n);
x = linspace(0,0.5,n); % create n linearly spaced elements bw 0 and 0.5

% iterate over the length of x (which in this case is n from above)
for i=1:length(x)
    y(i)=(Co/(2*D*t))*exp(-((x(i)^2)/(4*D*t)));
end

If you are unsure of the inputs for linspace or any other MATLAB command, just type (for example) help linspace in the Command Window for details/info on that command.

The trick now is to determine whether what is plotted is correct. I guess it all depends on the input values for the time and the concentration...

Answer 1

[Geoff Hayes]

Mo - A couple of things I did notice in the above code: there will be an an error/failure when the code tries to add the labels to the x and y axes of your plots - the commands are xlabel and ylabel respectively.

Also, consider creating an x vector using linspace to size x to be the same as y (for values from 0 to 50 centimetres). Then iterate over the length of x and plug x(i) into the equation. With an actual vector x, you can use this in your plot command as plot(x,y). Please check its documentation (type help plot) to see why this is different from what you have.