Different sizing on LHS and RHS

Hi there, I'm wondering if anyone could explain why this error is happening. My code is not solving the full 100 iterations when solving for Ma1. It only reaches a value of i=98 in the for loop solving for ma1. Is there a reason it is stopping? or is there a way to make it solve all 100? I don't understand why it has an upper cap.

When I set TL to 1.8, the solution runs perfectly fine, but when I move it to 1.9 or higher it does not solve all the values of Ma1.

%set variables 
To = 330;
MF_max = 0.2;
L_max= 2.1;
G = 1.4;
R = 287;
f = 0.003;
Po = 4;
TL = 2; % we can change this one 
%PSW = 0.4; % location of pipe in realtion to total pipe length (m)
PSW = linspace(0.36,0.69);  % min distance of shockwave in realtion to total pipe length
AL = length(PSW);
%PSW = 0.69  % max distance of shockwave in realtion to total pipe length
L12 = PSW*TL;
D = 0.03;  % we can change this one too
Ma4 = 1;
P4 = 1.013;

Lc3 = TL-L12;
FLcD3 = 4*f*Lc3/D;

syms ma3
Ma3a = [];
for i = 1:1:AL
    eqn1 = (-1./G).*(1-ma3.^-2)+((G+1)./(2.*G)).*log((((G+1)./2).*ma3.^2)./(1+((G-1)./2).*ma3.^2)) == FLcD3(i);
    Ma3a(1,i) = vpasolve(eqn1,ma3,0.01);
end
Ma3 = Ma3a;

syms ma2
Ma2a = [];
for i = 1:1:AL
    eqn2 = ((1+((G-1)./2).*ma2.^2)./(G.*ma2.^2-((G-1)./2))).^0.5 == Ma3(i);
    Ma2a(1,i) = vpasolve(eqn2,ma2,[1 100]);
end
Ma2 = Ma2a;

FLcD2 = (-1./G).*(1-Ma2.^-2)+((G+1)./(2.*G)).*log((((G+1)./2).*Ma2.^2)/(1+((G-1)./2).*Ma2.^2));

FL12D = 4.*f.*L12./D;

FLcD1 = FLcD2 + FL12D;

syms ma1
Ma1a = [];
for i = 1:1:AL
eqn3 = (-1./G).*(1-ma1.^-2)+((G+1)./(2.*G)).*log((((G+1)./2).*ma1.^2)./(1+((G-1)./2).*ma1.^2)) == FLcD1(i);
    Ma1a(1,i) = vpasolve(eqn3,ma1,1.1);
end
Ma1 = Ma1a;

Answer 1

[Alan Stevens]

It works ok using fzero, as follows:

%set variables 
To = 330;
MF_max = 0.2;
L_max= 2.1;
G = 1.4;
R = 287;
f = 0.003;
Po = 4;
TL = 2; % we can change this one 
%PSW = 0.4; % location of pipe in realtion to total pipe length (m)
PSW = linspace(0.36,0.69);  % min distance of shockwave in realtion to total pipe length
AL = length(PSW);
%PSW = 0.69  % max distance of shockwave in realtion to total pipe length
L12 = PSW*TL;
D = 0.03;  % we can change this one too
Ma4 = 1;
P4 = 1.013;
FL12D = 4.*f.*L12./D;
Lc3 = TL-L12;
FLcD3 = 4*f*Lc3/D;
FLcD2 = zeros(size(FLcD3));
FLcD1 = zeros(size(FLcD3));

eqn1 = @(ma3,i)(-1./G).*(1-ma3.^-2)+((G+1)./(2.*G)).*log((((G+1)./2).*ma3.^2)./(1+((G-1)./2).*ma3.^2)) - FLcD3(i);
eqn2 = @(ma2,ma3)((1+((G-1)./2).*ma2.^2)./(G.*ma2.^2-((G-1)./2))).^0.5 - ma3;
eqn3 = @(ma1, FLcD1)(-1./G).*(1-ma1.^-2)+((G+1)./(2.*G)).*log((((G+1)./2).*ma1.^2)./(1+((G-1)./2).*ma1.^2)) - FLcD1;

ma3 = zeros(1,100);
ma2 = zeros(1,100);
ma1 = zeros(1,100);
ma3_0 = 0.5;
ma2_0 = 2;
ma1_0 = 0.5;

nbrs = 1:100;
for i = nbrs
    
     ma3(i) = fzero(@(ma3) eqn1(ma3,i), ma3_0); 
     ma2(i) = fzero(@(ma2) eqn2(ma2, ma3(i)), ma2_0);
     FLcD2(i) = (-1./G).*(1-ma2(i).^-2)+((G+1)./(2.*G)).*log((((G+1)./2).*ma2(i).^2)/(1+((G-1)./2).*ma2(i).^2));
     FLcD1(i) = FLcD2(i) + FL12D(i);
     ma1(i) = fzero(@(ma1) eqn3(ma1,FLcD1(i)),ma2_0);
     ma3_0 = ma3(i);
     ma2_0 = ma2(i);
     ma1_0 = ma1(i);
     
end

plot(nbrs,ma3,nbrs,ma2,nbrs,ma1),grid
legend('ma3','ma2','ma1')

[Samuel Casallas]

Thanks a lot, that fixed my issue. Could you explain what ma3_0 does?

[Alan Stevens]

fzero needs an initial guess for the parameter it is trying to find. ma3_0 is the initial guess for ma3 (similarly ma2_0 is the initial guess for ma2, etc.).

[Samuel Casallas]

Just to clarify, eqn1, whatever it is, must always be the variable that is solved with different values within a range, in this case FLcD3, with the supsequent equations using values from the previous equations?

[Alan Stevens]

The way you define your equations you have eqn1 representing ma3, and eqn3 representing ma1. So

for eqn1, ma3 depends on G and FLcD3

for eqn2, ma2 depends on G and ma3

for eqn3, ma1 depends on G and FLcD1

I've just noticed that my listing above has the following incorrect line

ma1(i) = fzero(@(ma1) eqn3(ma1,FLcD1(i)),ma2_0);

This should be

ma1(i) = fzero(@(ma1) eqn3(ma1,FLcD1(i)),ma1_0);

i.e. it should use the initial value for ma1_0

This won't make any practical difference in this case as the initial values are just guesses anyway. However, it is best to put this right!

[Samuel Casallas]

Thanks a lot for your help. It makes more sense now and the calculation is more efficient than I had originally written.