Is it possible to fill a vector with my iterative results?

My script is set up as follows:

C0 = 3;

C = C0;

D = (1.16*(10^-5)*0.01^2)/60;

w = 15*10^-6;

n = 17000;

R = 100*10^-6;

Q = 0.003;

x = 0.25;

pass = 0;

f =(C*exp((-2*D*pi*R*x*n)/(w*Q)))

deltaC = C - f

scatter(pass, C)

hold on

while C > C0*0.3

C = C - deltaC;

f = (C*exp((-2*D*pi*R*x*n)/(w*Q)));

deltaC = C - f;

pass = pass +1;

scatter(pass, C)

end

hold off

basically it currently plots each iteration on a scatter graph, but I would rather it filled out a vector to be plotted instead so I can run sensitivity analysis on each variable and compare plots on the same graph

Answer 1

[William Rose]

@Mathew McFarlane,

If you do algebra on the equations in your code, you can show that

C=C0*K^p;

where p=pass number, and

K=exp(-2*D*pi*R*x*n/(w*Q));

You can calculate all the values of C in a single step, without a for loop, once you know pmax, the total number of passes needed. This is not obvious, since the increment df changes with every pass. However, we recognize that, although the increment df changes, the fractional change is the same on every pass: C is multiplied by K each time. We start with C=C0 and we want to make passes as long as C>0.3*C0. Therefore the number of passes needed, pmax, satisfies

C0*K^pmax > 0.3*C0

Therefore we divide both sides by C0, then take the log of both sides, then solve for pmax, to get

pmax=floor(log(0.3)/log(K));

Now we can get rid of the while loop, and compute all the values of C with one statement:

C=C0*K.^(1:pmax);

The plot shows C versus pass number, where pass number p=1:pmax.

I also added a plot of C0*exp(pass*ln(K)), to show that it is equivalent. The equivalence of the two plots demonstrates that one can thnk of this program as a plot of the exponential function.

The full code is below:

%MatthewMCode2.m
D = (1.16*(10^-5)*0.01^2)/60;
w = 15*10^-6;
n = 17000;
R = 100*10^-6;
Q = 0.003;
x = 0.25;
K=exp(-2*D*pi*R*x*n/(w*Q));     %compute K using constants above
pmax = floor(log(0.3)/log(K));  %compute max passes required
C0 = 3;
C=C0*K.^(1:pmax);               %compute vector C in one step

%plot results
figure;
subplot(2,1,1);
scatter(1:pmax,C);              %plot pass number versus C
grid on;
xlabel('Pass number'); ylabel('C');
subplot(2,1,2);
plot(1:pmax,C0*exp([1:pmax]*log(K)));  %plot of an alternate way to compute C
grid on;
xlabel('Pass number'); ylabel('C=exp(pass*ln(K))');

[Mathew McFarlane]

This is awesome, thank you so much. the explanations really helped as well, because I am far from proficient with this software.

[William Rose]

You're welcome@Mathew McFarlane.

Answer 2

[Mathieu NOE]

hello

simply index your output variable (C) - I also tweaked a bit the code in some other areas...

clc
clearvars

C0 = 3;
D = (1.16*(10^-5)*0.01^2)/60;
w = 15*10^-6;
n = 17000;
R = 100*10^-6;
Q = 0.003;
x = 0.25;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
pass = 1;
C(pass) = C0;
new_constant = exp((-2*D*pi*R*x*n)/(w*Q)); % defined once as this complex product does not change during the main loop (so avoid repeat it)
f =C*new_constant;
deltaC = C - f;
while C > C0*0.3
    pass = pass +1;
    C(pass) = C(pass-1) - deltaC;
    f = C(pass)*new_constant;
    deltaC = C(pass) - f;
end
scatter(1:pass, C)

[William Rose]

@Mathew McFarlane and @Mathieu NOE

@Mathieu NOE astutely reocognized that new_constant should be computed before entering the while loop. This quantity is called K in my proposed answer. In the code of @Mathieu NOE, vector C changes size on every pass thorugh the while loop. This is not recommended, because Matlab has to copy the previous value of the entire vector each time. A preferred and faster approach is to preallocate the vector. To pre-allocate, you need to know how long the vector is going to be, i.e. how many passes there will be through the while loop. See my suggeted answer for that. Or you can delete the while loop, and compute C with a single statement, as in my suggested answer.

Answer 3

[William Rose]

@Mathew McFarlane,

You compute a constant inside the loop which is the same every time, so compute it outside the loop to save time. I have done this in the attached code, and I have alled the constant K.

Also, yur algorithm to compute C canbe simplified considerably. Your vaue C is simply multipled by K each time through the loop. K is a little biy less than 1 so C gets smaller by a constant percentage with each pass. IN otherwords C delcines expnenetially.

I have to run, so I have not yet put it into a vector. We wanto precompute the size of C and eliminate the For loop altogether. Need to run now.

See attached code wihcih is simpler and faster than yours but give same output.

%MatthewMCode2.m
C = 3;
D = (1.16*(10^-5)*0.01^2)/60;
w = 15*10^-6;
n = 17000;
R = 100*10^-6;
Q = 0.003;
x = 0.25;
K=exp(-2*D*pi*R*x*n/(w*Q));
pass = 0;
scatter(pass, C)
hold on
while C > C0*0.3
    C = K*C;
    pass = pass +1;
    scatter(pass, C)
end
hold off