Different sizing on LHS and RHS

Question

Samuel Casallas 2021-3-27

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/785271-different-sizing-on-lhs-and-rhs

评论： Samuel Casallas 2021-3-28

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;

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Alan Stevens 2021-3-27

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/785271-different-sizing-on-lhs-and-rhs#answer_660016

在 MATLAB Online 中打开

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')