Different sizing on LHS and RHS

3 Ansichten (letzte 30 Tage)
Samuel Casallas
Samuel Casallas am 27 Mär. 2021
Kommentiert: Samuel Casallas am 28 Mär. 2021
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;

Akzeptierte Antwort

Alan Stevens
Alan Stevens am 27 Mär. 2021
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')
  5 Kommentare
Alan Stevens
Alan Stevens am 28 Mär. 2021
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
Samuel Casallas am 28 Mär. 2021
Thanks a lot for your help. It makes more sense now and the calculation is more efficient than I had originally written.

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)

Community Treasure Hunt

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

Start Hunting!

Translated by