Finite Difference Iteration Code Issues

I'm taking a groundwater course and the prof assumes that we know how to use matlab. I do not. He teaches material related to groundwater in class, and then tells us to model it in matlab, which he does not teach. My inability to use matlab is going to seriously impact how much of the material I learn in the class. I purchased the student version but have not been having much luck learning on my own. It is a simple finite difference method that we have to run until conversion. I have solved the problem already in excel but doing all of my assignments in excel is not going to help me learn to use matlab.

I would really appreciate if someone could help me out. In addition to solving, we also have to make a plot of the solution for C versus y at specified x values and also a table of C versus x at specified y values. I haven't even attempted that part yet.

clear;
%Define Parameters
r=0.5;
dt=16*r/0.119;
n=7;
 t=60;
%Initial Value of Function C
for i=2:n-2
for j=2:100
c(i,j-1)=0.2
end
%Value at boundary
for j=1:100
c(1,j)=0.;
c(n-1,j)=0.1;
end
%implement explicit method
for j=1:100
for i=1:n
c(i,j+1)=c(i,j)+r*(c(i+1,j)-2*c(i,j)+c(i-1,j));
end
clast=c %save the last guess
err(c)=max(abs(c(:)-clast(:)));
if err(C)<0.001
break
end

[Walter Roberson]

You do not appear to have a "y" in your code to plot C against.

[Roger Stafford]

It is very difficult to know how to help you with your problem. The finite difference procedure you are carrying out in the "%implement explicit method" part looks vaguely like an approximation to a partial differential equation of the form dc/dx = r*d(dc/dy)/dy with given boundary conditions on the left edge. That is just a guess. However, there are a number of errors present and the code is very incomplete so we can't really understand the details.

I would suggest that there is no "royal road" to understanding the matlab language. You simply have to devote enough time to learn its fundamentals if you are to write reasonably good code in it. There is no shortcut to doing that. All of us have had to go through that process. If you understand excel you should be well on your way.

[H]

The question that needs solving is a differential equation for diffusion with a boundary condition at both the top and bottom and i need it to iterate over time until steady state is reached. It needs to have 5 nodes in the i direction (I previously referred to as x) and iterate with time in the j direction.

I usually learn best by example but I can't find any clear examples of code that do what I am trying to do that I can pick apart and figure out what each part does.

If you could give me specific examples of what is wrong/incomplete it would really help me in fixing it.

Basically the question is there is a 20 cm long column that is broken down into 4 cm increments (so 5 nodes). The top has a concentration of 0 as a BC, the bottom has a concentration of 0.1 as the BC, the initial value for the column is 0.2, and it needs to reach steady state. We were asked to replace the D*dt/(dx*dx) term with varying values to show the difference it makes in the result (I've attempted to replace it with r in the above code).

[H]

Also it needs to iterate until the value of c at each node i at time j is within 0.001 of the value at the previous j.

[Roger Stafford]

Here are some presumed errors or inefficiencies I see in your code. If I were to try a little harder I could probably spot a few more. The nature of these tells me that you need to do some further study of matlab language fundamentals.

1. Missing 'end' after

 for i=2:n-2
   for j=2:100
     c(i,j-1)=0.2
   end

2. Another possible missing 'end' after

for j=1:100
  for i=1:n
    c(i,j+1)=c(i,j)+r*(c(i+1,j)-2*c(i,j)+c(i-1,j));
  end

3. In the inner loop in 2. c(0,j), c(7,j), and c(8,~) are being accessed, but c(0,j) is invalid and rows c(7,:) and c(8,:) have not been defined. I believe it ought to be:

     for i = 2:n-2
       c(i,j+1)=c(i,j)+r*(c(i+1,j)-2*c(i,j)+c(i-1,j));
     end

4. In 1. you filled c(2:n-2,1:99) with .2 values, but all columns after the first one are overwritten. Not an error, but a useless activity.

5. At "clast=c %save the last guess" you put the entire 6-by-100 matrix c into 'clast'. I seriously doubt that was your intention. That would make c(:)-clast(:) all zeros.

6. I haven't been able to figure out for sure what you are trying to do with the

if err(C)<0.001
  break
end

code. Presumably you are testing to see if the c column has reached a nearly steady state, but you appear to be comparing it with itself, not the previous c column. Something important is missing in your code there.

7. For efficient computation you should initially allocate the total size of matrix c you plan to fill with the command:

 c = zeros(n-1,100);

Note that that would make the setting of zeros in the first row unnecessary. This is also not an error but it can lead to very slow operation as columns are added one-by-one to 'c'.

8. To fill the sixth row with .1 values, you can simply write

 c(6,:) = .1

after doing the allocation in 7.

9. You can replace

for i = 2:n-2
  c(i,j+1)=c(i,j)+r*(c(i+1,j)-2*c(i,j)+c(i-1,j));
end

with

c(2:n-2,j+1) = c(2:n-2,j)+r*(c(3:n-1,j)-2*c(2:n-2,j)+c(1:n-3,j));

[H]

Thank you so much I really appreciate it. I will spend some time working on this. I understand I need to spend a lot more time learning matlab. My intention was to learn as much as I could in one weekend before my groundwater class next week so that I don't miss out on learning stuff about groundwater in the class due to my non-existent matlab skills. Unfortunately it was not listed as a prerequisite for the course to have knowledge of matlab but the professor is teaching it with the assumption that we know how to use matlab. I don't want to fall behind in the class because the groundwater related material will be important in my job. Matlab not so much aside from the fact that a working knowlege of it will be necessary to get the most out of this class.

Sorry to be so long-winded I just wanted to explain why it probably appears that I'm being lazy looking for an easy way to learn about matlab in an unreasonably short period of time.

[H]

For this part:

clast=c %save the last guess
err(c)=max(abs(c(:)-clast(:)));
if err(C)<0.001
break

My intention was to stop the loop once it reached almost steady state - when the value in each row of a column was within 0.001 of the value in that row of the previous column. Does this make more sense?

err(c)=max(abs(c(2:n-2,j)-c(2:n-2,j-1)));
  if err(c)<0.001
  break
  end

[Roger Stafford]

Yes it makes sense but that is not what you achieved by writing "clast=c". The effect of that line is, as I have said, to copy the entire 'c' matrix with all its 100 columns over to 'clast'. I think what you want is:

 clast = c(:,j);
 cpresent = c(:,j+1);
 if max(abs(clast-cpresent))<0.001
  break
 end

[H]

Everything seems to be functioning for me now with the exception of the break. I expected that if the if statement was met at say column 32 it would only display the matrix with 32 columns rather than continue to 100. However, it continues to display results up to column 100.

Is that what is supposed to happen?