Output for while loop

5 views (last 30 days)
Mateusz Brzezinski
Mateusz Brzezinski on 29 Jan 2020
Hello,
I wrote a script that gives roots based on the Newton Raphson iteration. However, it is looped for changing one factor (in this case A0). I would like to store all outputs (uA and uB) and correlate it with changing values of A0 so in text file I will have sth. like this:
A0=
0.01
uA=
0.06
uB=
0.83
A0=
0.02
uA=
0.09
uB=
0.82
and so on...
How I should do this?
This is my simple script:
diary Logs
format long
A0Start=0.000;
A0End=1.000;
A0Step=0.01;
for h=A0Start:A0Step:A0End
A0=h;
prec=1e-10;
uA=0.2;
uB=0.8;
B0=0.550;
KA=64.54860696;
KB=7.717484273;
KAB=96.50371794;
x=[uA;uB];
dev=1;
indx=0;
while dev>prec
f=[(A0*x(1)^2)+(((A0/(2*KA))^(1/2))*x(1))+(KAB/2)*(((A0*B0)/(KA*KB))^(1/2))*x(1)*x(2)-A0; (B0*x(2)^2)+(((B0/(2*KB))^(1/2))*x(2))+(KAB/2)*(((A0*B0)/(KA*KB))^(1/2))*x(1)*x(2)-B0];
J=[2*A0*x(1)+(A0/(2*KA))^(1/2)+(KAB*x(2)*((A0*B0)/(KA*KB))^(1/2))/2, (KAB*x(1)*((A0*B0)/(KA*KB))^(1/2))/2; (KAB*x(2)*((A0*B0)/(KA*KB))^(1/2))/2, 2*B0*x(2)+(B0/(2*KB))^(1/2)+(KAB*x(1)*((A0*B0)/(KA*KB))^(1/2))/2];
delx=-inv(J)*f;
x=x+delx;
dev=max(abs(f));
indx=indx+1;
end
disp 'for A0';
disp (h)
%disp 'last additive';
%disp(delx)
%disp 'iterations';
%disp (indx)
disp 'uA and uB found as';
disp (x)
end
diary off

Accepted Answer

Geoff Hayes
Geoff Hayes on 29 Jan 2020
Mateuz - try storing the data in an array which you can then (for example) write to file after the loop has completed. For example,
format long
A0Start=0.000;
A0End=1.000;
A0Step=0.01;
myData = zeros(length(A0Start:A0Step:A0End),3);
k = 1;
for h=A0Start:A0Step:A0End
A0=h;
prec=1e-10;
uA=0.2;
uB=0.8;
B0=0.550;
KA=64.54860696;
KB=7.717484273;
KAB=96.50371794;
x=[uA;uB];
dev=1;
indx=0;
while dev>prec
f=[(A0*x(1)^2)+(((A0/(2*KA))^(1/2))*x(1))+(KAB/2)*(((A0*B0)/(KA*KB))^(1/2))*x(1)*x(2)-A0; (B0*x(2)^2)+(((B0/(2*KB))^(1/2))*x(2))+(KAB/2)*(((A0*B0)/(KA*KB))^(1/2))*x(1)*x(2)-B0];
J=[2*A0*x(1)+(A0/(2*KA))^(1/2)+(KAB*x(2)*((A0*B0)/(KA*KB))^(1/2))/2, (KAB*x(1)*((A0*B0)/(KA*KB))^(1/2))/2; (KAB*x(2)*((A0*B0)/(KA*KB))^(1/2))/2, 2*B0*x(2)+(B0/(2*KB))^(1/2)+(KAB*x(1)*((A0*B0)/(KA*KB))^(1/2))/2];
delx=-inv(J)*f;
x=x+delx;
dev=max(abs(f));
indx=indx+1;
end
myData(k,1) = A0; % or h as they are the same
myData(k,2:3) = x; % assuming that x is a 1x2 array
k = k + 1;
end
  2 Comments
Mateusz Brzezinski
Mateusz Brzezinski on 29 Jan 2020
Geoff thank you a lot I will try it in free time,
Best.
Mateusz Brzezinski
Mateusz Brzezinski on 31 Jan 2020
Works like a charm.

Sign in to comment.

More Answers (0)

Products


Release

R2019a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!