How to extract a for loop output after each iteration?

Dear experts,

This is a for loop I am trying to implement. I want to store the result for u after each iteration of k1, c1, and rho1. I cannot figure out a way to do it.

Can some help with that?

Thank you!

for j =1:length(k1) % k1 is a matrix 1x10

for M =1:length(c1) % c1 is a matrix 1x10

for p = 1:length(rho1) % rho1 is a matrix 1x10

%u will be a matrix 17x1000

for t =1:Nt-1

u(1,t+1)= u(1,t) + 2*k1(j)*dt*(u(2,t)-u(1,t))/(rho1(p)*c1(M)*dx1^2) + 2*(hc(t)*(Ta(t) -u(1,t))+q(t))*dt/(rho1(p)*c1(M)*dx1);

for i=2:Nx-1

u(i,t+1) = u(i,t) + (k1(j).*dt)/(rho1(p).*c1(M).*dx1^2)*(u(i+1,t) - 2*u(i,t) + u(i-1,t));

end

u(Nx,t+1) = u(Nx-1,t);

end

end

end

[Rik]

Your code is not formatted, so it is difficult to read. Please edit your post.

It is hard to follow your code, but it looks like every iteration is overwriting the last.

You clearly understand indexing, why don't you use a 3D cell array to store the results? Alternatively, you can use save. What did you try?

Also, code without comments and with extremely short variable names become utterly worthless in a day or so. Nobody else can understand your code, and if you stop working on it over the weekend, you won't be able to understand it either. Do yourself a favor and explain to yourself what is happening and why. (relatedly, the code you posted doesn't show Nx is 17)

[Dyuman Joshi]

You could vectorize the inner most loop.

Also, have you pre-allocated the variable u?

"I want to store the result for u after each iteration of k1, c1, and rho1."

That would be a time consuming task.

Are you sure you want to save the values after each iteration and not after all the iterations?

What is the objective here? It would be helpful if you can specify what you are trying to do and provide relevant details.

Answer 1

[Star Strider]

Perhaps —

u(j,m,p,t) = ... ;

in an appropriate place in the loop.

This matrix could be huge and difficult to deal with.

Consider using the squeeze function if you want to easily eliminate any singletion dimensions that may arise if you are processing it later.

.

[Sanley Guerrier]

Yes, you are right Star Strider. It returns a huge matrix, I cannot deal with.

[Star Strider]

The problem may be the length of ‘t’. Consider carefully how long it actually needs to be to do what you want. Probably shorter is better here.

[Sanley Guerrier]

The length of t is very long and that is how I want it. I thought there were a simplest way to do it without implementing conditional logic.

[Rik]

If the results don't fit in a single array, you will have to save the intermediate steps in a mat file.

What do you intend to do next? That might help us help you to figure out if a more effective method exists.

[Sanley Guerrier]

The objective is to conduct a sensitivity analysis for each iteration and then compare them. Maybe I didn't use the right approach.

[Star Strider]

A sensitivity analysis (at least in my experience) is defined as calculating the partial derivative of a function (usually an objective function) with respect to each parameter. You can do this numerically (using the gradient function) or symbolically using the Symbolic Math Toolbox, I would go with the symbolic approach, however that is simply how I would deal with this sort of problem.

[Sanley Guerrier]

I already solved the partial derivative numerically by guessing k1, c1, rho1. Now I want to try different combination of k1, c1, and rho1 and find the best scenario that will probably optimize the function.

[Star Strider]

I would not suggest guessing at the values. Solve them symboilcally instead. (The Symbolic Math Toolbox should make that straightforward.) Unless the function has a step of some sort (so the derivative at that point is infinite), it should be continuously differentiable.

The gradient is likely provided by the Jacobian function of the objective function (again, symbolic). That is used to calculate the gradient descent, resulting the the most optimal parameter set.

From the symbolic calculations, you can create functions that can be evaluated numerically by using the matlabFunction function. This should make everything much easier.

[Sanley Guerrier]

Thank you for the suggestion. Solve them symbolically sounds great, I will try that.

[Star Strider]

As always, my pleasure!

[Sanley Guerrier]

Hi Star Strider,

I followed your advice, I try to optimize the parameters k, c, and rho using fminunc. But there is something wrong with the ojective function or the fminunc.

Could you expertise my code, please?

%Objective function

function obj = globalfun(x)

k =x(1);

c= x(2);

rho =x(3);

u= compute_Temp(k, c, rho); %Model depending on parameters k, c, and rho

data = readtable("data.xlsx");

T = data{:, 2:17};

difference = T(:,2:end)' - u(2:end,:);

obj = sqrt( sum(difference(:).^2)/(length(difference(:))-1) );

end

%% unconstrained optimization file

%Parameter initial guess

k = 1.2;

c = 900/84600;

rho = 2360;

x0 = [k; c; rho];

fun = globalfun; %function handle

[x, fval] = fminunc(fun, x0);

%optimization parameter values

k = x(1);

c = x(2);

rho = x(3);

disp(['k: ' num2str(k)])

disp(['c: ' num2str(c)])

disp(['rho: ' num2str(rho)])

%calculate model with update parameter

x = [k; c; who];

u = compute_Temp(k, c, rho);

[Star Strider]

There are several problems.

First, the ‘compute_Temp’ function is nowhere to be seen. Neither is ‘data.xlsx’ however unless the Run feature is back up and running (it wasn’t all day yesterday and still isn’t when I just now checked), I can’t run your code here anyway, so we need to wait for all those to appear (or reappear) before I can do anything.

Second, the readtable call belongs in the script, not the objective function, because it will be incredibly slow. Pass whatever values you need from it to the objective function as an extra parameter.

Third, there is no need to do detailed calculations inside the objective function. It only needs to return the norm of ‘difference’ and since ‘compute_Temp’ seems to be the actual objective function anyway, the fminunc call could just use it (except for the restriction of only using the first 17 rows or columns of it).

A more efficinet version coould be:

data = readtable("data.xlsx");

T = data{:, 2:17};

k = 1.2;

c = 900/84600;

rho = 2360;

x0 = [k; c; rho];

fun = @(x) globalfun(x,T); %function handle

[x, fval] = fminunc(fun, x0);

%optimization parameter values

k = x(1);

c = x(2);

rho = x(3);

% disp(['k: ' num2str(k)])

% disp(['c: ' num2str(c)])

% disp(['rho: ' num2str(rho)])

%calculate model with update parameter

% x = [k; c; who];

fprintf('\nk = %8.4f\nc = %8.4f\nrho = %8.4f\n',k,c,rho)

u = compute_Temp(k, c, rho);

function obj = globalfun(x,T)

k =x(1);

c= x(2);

rho =x(3);

u= compute_Temp(k, c, rho); %Model depending on parameters k, c, and rho

obj = norm(T(:,2:end)' - u(2:end,:));

end

Try this.

.

[Sanley Guerrier]

It works fine now. Thank you very much!

why do you replace the:

obj = sqrt( sum(difference(:).^2)/(length(difference(:))-1) );

by

obj = norm(T(:,2:end)' - u(2:end,:));

what is the main difference?

[Star Strider]

As always, my pleasure!

I replaced it because you only need the norm of ‘difference’. It’s more efficient, and therefore faster. It doesn’t need to be scaled or normalised at that point in your code.

.

[Sanley Guerrier]

Yes, you're right. I tried both and the result clearly shows that the norm of difference is more efficient and faster.

[Star Strider]

Thank you!