Why am I getting an imaginary number in my loop?

Hi everyone,

For my loop I keep getting imaginary numbers. This is not supposed to happen and was wondering if anyone had any tips. In the code, theta_o is supposed to converge when it equals 0.1696, but when it goes negative that makes the convergence not happen

%% variables
Ctreq = 0.008;
Ddo = 0.01;
Clalpha = 5.76;
sigma = 0.1;
Nb = 4;
dr = 0.025;

%% initalize
Ct = 0;
Cpi = 0;
theta_tw = 0;
theta_j = 0.1;


theta_o = 6.*Ctreq./sigma./Clalpha - 3./4.*theta_tw + 1.0607.*sqrt(Ctreq);
theta_o = [theta_o];
    
        for r = 0:0.025:1
            theta_o = theta_o(end);
            lambda = (sigma.*Clalpha./16).*(sqrt(1+32.*theta_o.*r./sigma./Clalpha)-1);
            
            dCt = (sigma.*Clalpha./2).*(theta_o.*r.^2 - lambda.*r).*dr;
            Ct = Ct + dCt;
            Cpi = Cpi + lambda.*dCt;
            
            theta_j = theta_o + (6.*(Ctreq-Ct)./sigma./Clalpha + 1.0607 .* (sqrt(Ctreq)-sqrt(Ct)));
            theta_o = [theta_o theta_j];
            
            if abs((theta_j-theta_o(1))/theta_j) > 0.001
                break 
            end 
                
        end 

Answer 1

[Walter Roberson]

When your dCt is positive, your Ct increases. If it does that often enough, it can end up being larger than Ctreq, and then (Ctreq-Ct) would be negatie and sqrt(Ctreq)-sqrt(Ct) would be negative, so theta_j would start decreasing. If your dCt does not start going negative then theta_j will get smaller and smaller until it becomes negative, and then the part inside your sqrt() for lambda is at risk of going complex.

%% variables
Ctreq = 0.008;
Ddo = 0.01;
Clalpha = 5.76;
sigma = 0.1;
Nb = 4;
dr = 0.025;

%% initalize
Ct = 0;
Cpi = 0;
theta_tw = 0;
theta_j = 0.1;


theta_o = 6.*Ctreq./sigma./Clalpha - 3./4.*theta_tw + 1.0607.*sqrt(Ctreq);
theta_o = [theta_o];
    
        for r = 0:0.025:1
            theta_o = theta_o(end);
            part1 = 1+32.*theta_o.*r./sigma./Clalpha
            lambda = (sigma.*Clalpha./16).*(sqrt(part1)-1);
            
            dCt = (sigma.*Clalpha./2).*(theta_o.*r.^2 - lambda.*r).*dr;
            Ct = Ct + dCt;
            Cpi = Cpi + lambda.*dCt;
            
            if Ctreq < Ct
              fprintf('Caution, watch out for negative, Ctreq = %g, Ct = %g\n', Ctreq, Ct);
            end
            theta_j = theta_o + (6.*(Ctreq-Ct)./sigma./Clalpha + 1.0607 .* (sqrt(Ctreq)-sqrt(Ct)))
            theta_o = [theta_o theta_j];
            
            if abs((theta_j-theta_o(1))/theta_j) < 0.001
                break 
            end 
                
        end 

part1 = 1

theta_j = 0.3564

part1 = 1.4950

theta_j = 0.5342

part1 = 2.4839

theta_j = 0.7108

part1 = 3.9615

theta_j = 0.8852

part1 = 5.9177

theta_j = 1.0564

part1 = 8.3362

theta_j = 1.2232

part1 = 11.1935

theta_j = 1.3842

part1 = 14.4574

theta_j = 1.5377

part1 = 18.0858

theta_j = 1.6819

part1 = 22.0242

theta_j = 1.8148

part1 = 26.2049

theta_j = 1.9338

part1 = 30.5441

theta_j = 2.0365

part1 = 34.9410

theta_j = 2.1199

part1 = 39.2766

theta_j = 2.1812

part1 = 43.4126

theta_j = 2.2173

part1 = 47.1928

theta_j = 2.2250

part1 = 50.4437

Caution, watch out for negative, Ctreq = 0.008, Ct = 0.00946112

theta_j = 2.2014

part1 = 52.9786

Caution, watch out for negative, Ctreq = 0.008, Ct = 0.0116324

theta_j = 2.1441

part1 = 54.6020

Caution, watch out for negative, Ctreq = 0.008, Ct = 0.0140133

theta_j = 2.0507

part1 = 55.1170

Caution, watch out for negative, Ctreq = 0.008, Ct = 0.0165538

theta_j = 1.9200

part1 = 54.3347

Caution, watch out for negative, Ctreq = 0.008, Ct = 0.0191842

theta_j = 1.7515

part1 = 52.0856

Caution, watch out for negative, Ctreq = 0.008, Ct = 0.021814

theta_j = 1.5458

part1 = 48.2334

Caution, watch out for negative, Ctreq = 0.008, Ct = 0.0243333

theta_j = 1.3051

part1 = 42.6905

Caution, watch out for negative, Ctreq = 0.008, Ct = 0.0266153

theta_j = 1.0330

part1 = 35.4337

Caution, watch out for negative, Ctreq = 0.008, Ct = 0.0285226

theta_j = 0.7350

part1 = 26.5197

Caution, watch out for negative, Ctreq = 0.008, Ct = 0.0299175

theta_j = 0.4181

part1 = 16.0969

Caution, watch out for negative, Ctreq = 0.008, Ct = 0.0306818

theta_j = 0.0909

part1 = 4.4079

Caution, watch out for negative, Ctreq = 0.008, Ct = 0.0307875

theta_j = -0.2377

part1 = -8.2453

Caution, watch out for negative, Ctreq = 0.008, Ct = 0.0301302

theta_j = -0.5575 + 0.0070i

part1 = -21.4553 + 0.2827i

Caution, watch out for negative, Ctreq = 0.008, Ct = 0.0282025

theta_j = -0.8513 + 0.0255i

part1 = -34.4693 + 1.0644i

Caution, watch out for negative, Ctreq = 0.008, Ct = 0.0249317

theta_j = -1.1004 + 0.0586i

part1 = -46.3800 + 2.5247i

Caution, watch out for negative, Ctreq = 0.008, Ct = 0.0203365

theta_j = -1.2859 + 0.1083i

part1 = -56.1509 + 4.8141i

Caution, watch out for negative, Ctreq = 0.008, Ct = 0.0145519

theta_j = -1.3888 + 0.1758i

part1 = -62.6515 + 8.0591i

theta_j = -1.3912 + 0.2629i

part1 = -64.6972 + 12.4165i

theta_j = -1.2805 + 0.3776i

part1 = -61.2485 + 18.3566i

theta_j = -1.0693 + 0.5269i

part1 = -52.4671 + 26.3428i

theta_j = -0.7767 + 0.6924i

part1 = -38.9156 + 35.5806i

theta_j = -0.4170 + 0.8511i

part1 = -21.0064 + 44.9207i

theta_j = -0.0248 + 0.6695i

part1 = -0.3417 + 36.2645i

theta_j = 0.3657 + 0.4467i

part1 = 21.3157 + 24.8181i

theta_j = 0.7302 + 0.2008i