For Loop Indices Manipulation

Hi all. I have a function file which has a for loop which, depending on the conditions, returns certain equations that should end up creating a vector of differential equations, of number x that is determined by user input, that would be solved with 'ode45' in the main m-file. The issue is that the loop must start at 3 in order to have real positive integers.

function [ DeltaO2Concentration ] = O2Solver_beta(t,O2Conc,c,O2Leak_term,nodes_num,slt_nodes2,last_node)
%UNTITLED Test function used to solve 'first attempt' code.
%   System of ODE diffusion equations that will be solved using 'ode45'
%   function. Inputs are constant composed of universal and calculated
%   constants, time (independent variable), O2Conc (dependent variable), 
%   O2Leak which is a separate term or something.

%VECTOR LENGTH WILL NEED TO VARY ACCORDING TO USER INPUT%

% DeltaO2Concentration = [0;
%  c*t*(O2Conc(1)-O2Conc(2)+O2Conc(3)-O2Conc(2));
%     c*t*(O2Conc(2)-O2Conc(3)+O2Conc(4)-O2Conc(3));
%     c*t*(O2Conc(3)-O2Conc(4)+O2Conc(5)-O2Conc(4))%+O2Leak;
%     c*t*(O2Conc(4)-O2Conc(5)+O2Conc(6)-O2Conc(5));
%     c*t*(O2Conc(5)-O2Conc(6)+O2Conc(7)-O2Conc(6));
%     c*t*(O2Conc(6)-O2Conc(7)+O2Conc(8)-O2Conc(7))+O2Leak;
%     c*t*(O2Conc(7)-O2Conc(8)+O2Conc(9)-O2Conc(8))+O2Leak;
%     c*t*(O2Conc(8)-O2Conc(9)+O2Conc(10)-O2Conc(9))+O2Leak;
%     c*t*(O2Conc(9)-O2Conc(10)+O2Conc(11)-O2Conc(10));
%     c*t*(O2Conc(10)-O2Conc(11)+O2Conc(12)-O2Conc(11));
%     c*t*(O2Conc(11)-O2Conc(12)+O2Conc(13)-O2Conc(12));
%     c*t*(O2Conc(12)-O2Conc(13)+O2Conc(14)-O2Conc(13))%+O2Leak;
%     c*t*(O2Conc(13)-O2Conc(14)+O2Conc(15)-O2Conc(14));
%     c*t*(O2Conc(14)-O2Conc(15)+O2Conc(16)-O2Conc(15))%+O2Leak;
%     c*t*(O2Conc(15)-O2Conc(16)+O2Conc(17)-O2Conc(16));
%     c*t*(O2Conc(16)-O2Conc(17)+O2Conc(18)-O2Conc(17));
%     c*t*(O2Conc(17)-O2Conc(18)+O2Conc(19)-O2Conc(18))%+O2Leak; 
%     c*t*(O2Conc(18)-O2Conc(19)+O2Conc(20)-O2Conc(19));
%     0];  

terminate= nodes_num;
DeltaO2Concentration= zeros(1,terminate);      
cnt= 1;
for i= 3:length(DeltaO2Concentration)
    if i==3 && ismember(i,slt_nodes2)
        DeltaO2Concentration(cnt) = c*t*(O2Conc(i-2)-O2Conc(i-1))+O2Leak_term;  %first node, leakage
    elseif i==3
        DeltaO2Concentration(cnt) = c*t*(O2Conc(i-2)-O2Conc(i-1));  %first node, NO leakage
    elseif i==last_node && ismember(i,slt_nodes2)
        DeltaO2Concentration(cnt) = c*t*(O2Conc(i)-O2Conc(i-1))+O2Leak_term;  %last node, leakage
    elseif i==last_node
        DeltaO2Concentration(cnt) = c*t*(O2Conc(i)-O2Conc(i-1));  %last node, NO leakage
    elseif ismember(i,slt_nodes2)
        DeltaO2Concentration(cnt) = c*t*(O2Conc(i-2)-O2Conc(i-1)+O2Conc(i)-O2Conc(i-1))+O2Leak_term;  %other nodes, leakage
    else
        DeltaO2Concentration(cnt) = c*t*(O2Conc(i-2)-O2Conc(i-1)+O2Conc(i)-O2Conc(i-1));  %other nodes, NO leakage
    end
    cnt= cnt+1;
end
DeltaO2Concentration= DeltaO2Concentration.';   %column vector for ODE45

nodes_num: number of nodes user inputs, which translates to # of ODE equations.

slt_nodes2: user input selection of specific nodes that will include the 'O2Leak_term' in equation.

The commented out stuff in the beginning is just a general template that I use to create the loop.

If I run this as it is, the 'DeltaO2Concentration' vector is the correct length, but I feel like Matlab just puts in zero vectors or something at the end because I think it should return a vector two less than the desired, but I'm not sure how to interpret it or fix it to where I get the correct output. Any help in appreciated.

[Adam]

terminate= nodes_num;
DeltaO2Concentration= zeros(1,terminate);

is what tells Matlab exactly what to put in the vector initially and that will remain in those elements you do not overwrite. If you simply want to have two fewer elements then just make the vector shorter by 2 and update your indexing accordingly.

for i= 3:( length(DeltaO2Concentration) + 2 )

would be one of numerous ways to achieve this.

[Areeb Hussain]

Thanks Adam. That was one of the ways I tried to fix the problem. However, in doing that, I also needed to 'add 2' to other variables that are impacted or impact the function file. And when all that was accomplished, 'Deltao2Concentration' was length two greater than the desired output. So I'm trying to figure out another way to go about fixing this.

[Adam]

Well, I guess you could just truncate the vector afterwards if that is easier. Since you are indexing into it with something that should reach 2 less than the size of the vector when the loop terminates it should never write to the last two elements of your vector so just

DeltaO2Concentration( end-1:end ) = [];

should do the trick equally.

I would probably index into that vector using i-2 personally, or at the very least set

cnt = i-2;

in the loop rather than incrementing cnt, but it should all amount to the same thing, just preference

[dpb]

Just how many ODEs do you have? The derivative vector must be of that length.

We can certainly clean up the code as written significantly, but in reading your description I'm just not following the actual problem definition to know which is the right solution.

[Steven Lord]

Did the technique I described my answer in one of your previous questions not work?

[Areeb Hussain]

Thanks Adam. That could possibly be a way to do it, but I don't want to delete any correct elements of the vector (it's hard for me to see which ones are and aren't correct because these get fed into an ode solver function. But I'll try to look into it.

[Adam]

Well, I mean trim them before they get fed in, straight after your for loop.

[Areeb Hussain]

dpb, the number of ODEs in the vector system depends on user input (# of nodes). If they input 20, the function file should return a system of 20 ODEs that are correctly configured. All the commented stuff was just me messing with it to see how to fix it.

[Areeb Hussain]

Steven, I really didn't know how to proceed with that option, so I figured a for loop would be easier, until I ran into this issue.

[dpb]

OK, so now we know you do want N ODEs for the input number.

What you've got builds 20 but two are empty which may/may not be what the actual ODE system needs.

What is RHS for the first entries supposed to be?

If you had a smaller system of say five, show us the actual ODEs to be solved to really nail down how to write this correctly for the system.

[Areeb Hussain]

So if the user wants 5 nodes, and they select nodes 2 and 3 to include the 'O2leak_term,' then the function should return:

DeltaO2Concentration= [c*t*(O2Conc(1)-O2Conc(2));
c*t*(O2Conc(1)-O2Conc(2)+O2Conc(3)-O2Conc(2))+O2Leak_term;
c*t*(O2Conc(2)-O2Conc(3)+O2Conc(4)-O2Conc(3))+O2Leak_term;
c*t*(O2Conc(3)-O2Conc(4)+O2Conc(5)-O2Conc(4));
c*t*(O2Conc(5)-O2Conc(4))]

But these equations can vary based on the number of nodes, and the user selected nodes for 'O2leak_term'.

[Areeb Hussain]

And Adam, yeah I knew what you meant. I tried that, but that gave me another issue, because when I use 'ode45' in the main file, I input an initial condition vector y0 the same size as the ODE vector, but the truncation somehow messes up this up when running 'ode45', saying the dimensions of DeltaO2Concentration and y0 need to be the same. I'm trying to fix that right now too.

[Guillaume]

Coming rather late in the conversation and not having much to say to answer the actual problem, I'll just comment on this:

terminate= nodes_num;
DeltaO2Concentration= zeros(1,terminate);      
for i= 3:length(DeltaO2Concentration)

What's wrong with:

DeltaO2Concentration= zeros(1, nodes_num);
for i = 3:nodes_num

Why rename the variable and make the reader do extra mental work in order to understand the bounds of the loop?

Answer 1

[dpb]

OK, from that definition, a looping solution.

function dO2Conc = O2Solver(t,O2Conc,c,O2Leak,nodes,slt_nodes)
  dO2Conc= zeros(nodes,1);      
  dO2Conc(1) = O2Conc(1)-O2Conc(2);  % first node
  for i= 2:nodes-1
    dO2Conc(i) = O2Conc(i-1)-2*O2Conc(i)+O2Conc(i+1);
  end
  dO2Conc(nodes) = O2Conc(nodes)-O2Conc(nodes-1);
  dO2Conc=c*t*dO2Conc;
  dO2Conc(slt_nodes)=dO2Conc(slt_nodes)+O2Leak;

It looks like the last node breaks the pattern, though??? Is that really right?

The above will, however, provide you with the required length of the RHS vector based on requested node number with no artificial zero entries as had before. The last two lines add the normalization and then apply the constant leakage term on the selected nodes via vector addressing. That avoids the special-casing in the main build.

I see where Steven's headed with the spdiag route but will have to scratch my head a little more over the linear algebra to directly apply it -- the middle terms are easy, incorporating the boundary conditions means modifying the example some and I don't have time to mess with it just now, sorry.

However, to leave some bread crumbs, write out and study the pattern of the subscripts versus row index for each of the terms and then compare to the indices returned in the example Steven pointed you at. Note the middle elements are all [1 -2 1] for three columns. See any parallel to your expansion above?

Oh, wait a minute--

>> n=5;
>> A = spdiags([e [-1; -2*e(2:end-1);-1] e], -1:1, n, n);
>> spdiags(A)
ans =
   1    -1     0
   1    -2     1
   1    -2     1
   1    -2     1
   0    -1     1
>>

Hmmm...almost, last row should be [-1 1 0] it the form as written above is correct. Which isn't symmetric with first node that made me ask if that is correct.

Gotta' run, though, leave as "exercise for the student". :)

[Areeb Hussain]

Thanks dpb, this solution is definitely better and simpler. It outputs what I wanted and works with the 'ode45' function. Thank you. Also, I now see where you were going with the 'spdiags' function, Steven. For some reason I didn't put two and two together. Thanks for the help guys. I might need some more in the future haha.

[dpb]

You're welcome; just needed to know the form precisely.

I see there's an error above in the application of the leak term, however; in copying the previous line I forgot to remove the c*t* scaling factor so it'll be applied twice for those nodes; remove the multiplier there; it's already been applied universally--I did edit the Answer to correct.