While loop within for loop problem

Hello all,

I am writing a simple code to calculate the electric field strength in the insulation of a cable directly from the Maxwell equations at each time step. I want to iterate every time step until error is less than 0.1% to make my solution accurate. The code is not working properly because the while loop iterates only once in the first iteration of the r loop and I have no idea why ... ?

Here is the part of code:

z = size(Y_temp,1);

for r = 2:z-1
    V(r,1) = 450 + r; V(r,end) = 0;
    while error3 > tolerance %&& error4 > tolerance && error5 > tolerance && error6 > tolerance
        V_old = V;
        for c = 2:numel(r1)-1       
            ch(r,c) = ch(r-1,c) - (t/(2*h1)) * (J(r,c+1) - J(r,c-1))    %charge from current continuity
            V(r,c) = (h1/(4*r1(i))) * (V(r,c+1) - V(r,c-1)) + (V(r,c+1) + V(r,c-1))/2 + ((h1^2)/(2*epsilon)) * ch(r,c);     %electric potential for Poissons equation
            E(r,c) = (-V(r,c+1) + V(r,c-1))/(2*h1);     %Irrotational electric field
            sigma(r,c) = sigma_0 * exp(alpha_T * Y_temp(r,c)) * exp(gamma_E * E(r,c));  %empiric equation of conductivity
            J(r,c) = sigma(r,c) * E(r,c);     %Ohm's law           
        end
        error3 = max(abs((V - V_old) ./V_old));
    end
end

Unrecognized function or variable 'z'.

Thank you very much in advance for your help.

[Torsten]

Each V on the right-hand side of the equations within the for-loop has to be replaced by V_old.

Maybe the same is true for some other iteration variables - I don't understand the equation(s) you are iterating.

Answer 1

[Swati Khese]

I think there is problem in this line "z = size(Y_temp,1);"

Maybe Y_temp is a vector with two rows, so z becomes 2 and for loop runs for only once.

Can you share size of Y_temp variable? [size(Y_temp)]

[Kamil]

Hi Swati,

You are right about that. Y_temp is a matrix of 86374 x 95 of the temperatures that were calculated earlier in the other Matlab script.

Answer 2

[Kamil]

The answer to the problem I was facing was quite simple. For each row iteration, the value of error needs to be updated. The loop only runs in the first iteration. This is because after the first iteration the value of error was less than < tolerance. I rearranged the code a little. I've done a vectorisation in column direction and I think I got an accurate solution, but unfortunately it takes ridiculus ammount of time. Matlab profiler shows that the 152609 s .... I'm wondering: Is there a way to dramatically speed up this code? here is a rearranged code:

clc;
%format long
% ----- Calculation of electric field strenght and conductivity ------------------------ %
% Input data
U = 450;                                             % Voltage kV
sigma_0 = (1e-16) / 1000;                            % conductivity at 0 celsuic deegre 
epsilon_r = 3.5 / 1000;                              % relative permittivity
epsilon_0 = (8.97e-12) / 1000;                       % dielectric permittivity of the free space
epsilon = epsilon_r * epsilon_0;                     % absolute permitivity
alpha_T = 0.1;%0.084;                                % temperature dependency coefficient [1/Celsius]
gamma_E = 0.03;%0.0645;                              % field dependency coefficient
r_inner = 26.5;                                      % inner insulation radius
r_outer = 47.5;                                      % outer insulation radius
nx1 = 30;
t = 1;
r1 = linspace(r_inner,r_outer,nx1);                  % preparing space for field calculation
radius = [0 r1 0];
tolerance = 1e-3;
error = inf;
h1 = (r_outer - r_inner)/nx1;
% ---- Calculating of initial conditions ---- %
V0 = zeros(1,nx1+2);
ch0 = zeros(1,nx1+2);
E0 = zeros(1,nx1+2);
J0 = zeros(1,nx1+2);
sigma0 = zeros(1,nx1+2);
% ----- initial conditions -----%
while error > tolerance | error2 > tolerance
    V0_old = V0;
    ch0_old = ch0;
    for i=2:numel(V0)-1
        E0(i) = (U)/ (radius(i) * (log(radius(end-1)/radius(2))));
        sigma0(i) = sigma_0 * exp(alpha_T * (Y_temp(1,i-1) - 273.15) + gamma_E * E0(i));
        J0(i) = sigma0(i) * E0(i);
    end
    for i = 2: numel(V0)-1
        V0(1) = U; V0(end) = 0; 
        %V0(i) = (V0(i+1) - V0(i-1))* (h1/(4*radius(i))) + (V0(i+1) + V0(i-1))/2;
        V0(i) = V0(i+1) * (h1/(4*radius(i)) + 0.5) + V0(i-1) * (0.5 - h1/(4*radius(i)));
    end
    for i=2:numel(V0)-1
        %ch0(i) = ((epsilon/(h1)) * (E0(i+1) - E0(i)));
        ch0(i) = (-epsilon/(2*h1)) * (V0(i+1) - V0(i-1)) + (-radius(i) * epsilon/h1^2) * (V0(i+1) + V0(i-1) - 2*V0(i));
        
    end
    error = max(abs((V0 - V0_old)));
    error2 = max(abs((ch0 - ch0_old)));
end

Unrecognized function or variable 'Y_temp'.

% ----- solution -----%
z = size(Y_temp,1);
E = zeros(z,nx1+2);
ch = zeros(z,nx1+2);
sigma = zeros(z,nx1+2);
J = zeros(z,nx1+2);
V = zeros(z,nx1+2);
t_relax = zeros(z,nx1+2);
E(1,:) = E0;
ch(1,:) = ch0;
sigma(1,:) = sigma0;
J(1,:) = J0;
V(1,:) = V0;
E_old(1,:) = zeros(1,numel(radius));
n = 1;

for r = 2:z
    error3 = inf;
    V(r,1) = 450; V(r,end) = 0;
    while error3 > tolerance   
        E_old = E;
        V_old = V;
        
        % ---- current continuity equation ----%
        ch(r,2) = ch(r-1,2) - (t/(2*h1)) .* (J(r,3) + J(r-1,3) - J(r,2) - J(r-1,2));
        ch(r,3:end-2) = ch(r-1,3:end-2) - (t/(4*h1)) .* (J(r,4:end-1) + J(r-1,4:end-1) - J(r,2:end-3) - J(r-1,2:end-3));
        ch(r,end-1) = ch(r-1,end-1) - (t/(2*h1)) .* (J(r,end-1) + J(r-1,end-1) - J(r,end-2) - J(r-1,end-2));
        % ---- Poisson's equation ----%
        V(r,2:end-1) =(-ch(r,2:end-1)./epsilon + V(r,1:end-2)./h1 - radius(2:end-1)./h1^2 .* (V(r,3:end) + V(r,1:end-2))) ./ (1/h1 - (2*radius(2:end-1))./h1^2);
        % ---- irrotatinal Electric field E = -div V ----%
        E(r,2:end-1) = (V(r,1:end-2) - V(r,2:end-1)) ./ (h1);
        % ---- Conductivity gradient ----%
        sigma(r,2:end-1) = (sigma_0 * exp(alpha_T .* (Y_temp(r,1:nx1) - 273.15) + gamma_E .* E(r,2:end-1)));
        % ---- Ohm's law ----%
        J(r,2:end-1) = sigma(r,2:end-1) .* E(r,2:end-1);
        
        error3 = max(abs((E(r,2:end-1) - E_old(r,2:end-1))));
        
        n = n+1;
    end
end
n

[Torsten]

I'm wondering: Is there a way to dramatically speed up this code?

Instead of using fixed-point iteration, use a standard nonlinear solver to solve for the unknowns. It reduces computation time drastically in most cases.

[Kamil]

Thanks a lot for your answer. Since I'm a coding novice, do you mean using fsolve for example? Or is there any other solver that you would recommend me to use?

[Torsten]

Yes, "fsolve" would be my first choice.

[Kamil]

fsolve is part of the optimisation toolbox. unfortunately I do not have any toolboxes, I only have clean Matlab. Are there other ways to do optimisation?

[Torsten]

Here is a very primitive Newton method. Or use "Octave".

f = @(x)[x(1)^2-3;x(1)*x(2)-3.5];
u0 = [1;1];
sol = nls(f,u0)

iter = 5

sol = 2×1

1.7321 2.0207

function u = nls(f,uold)
  u = zeros(numel(uold),1);
  itermax = 30;
  eps = 1e-6;
  error = 1.0e5;
  uiter = uold;
  iter = 0;
  while error > eps && iter < itermax
    u = uiter - Jac(f,uiter)\f(uiter);
    error = norm(u-uiter);
    iter = iter + 1;
    uiter = u;
  end
  iter
end

function J = Jac(f,u)
  y = f(u);
  h = 1e-6;
  for i = 1:numel(u)
    u(i) = u(i) + h;
    yh = f(u);
    J(:,i) = (yh-y)/h;
    u(i) = u(i) - h;
  end
end