problem in the for - loop

2 次查看(过去 30 天)
Gleb
Gleb 2021-1-20
评论: Gleb 2021-2-7
Here is the problem
Ive got a FEM code
function main1
clc
n = 100;
x0 = ones(3*n,1);
sol = fsolve(@(x)fun(x,n),x0);
norm(fun(sol,n))
x = ((1:n)-1)/(n-1);
plot(x,sol(1:n))
end
function res = fun(z,n)
eta=1.0;
beta = 1.0;
x = ((1:n)-1)/(n-1);
dx = 1/(n-1);
y1 = z(1:n);
y2 = z(n+1:2*n);
y3 = z(2*n+1:3*n);
for i=1:length(y1)
Y1=y1(i);
end
for i=1:length(y2)
Y2=y2(i);
end
res_y1 = zeros(n,1);
res_y2 = zeros(n,1);
res_y3 = zeros(n,1);
res_y1(1) = y1(1)-10.0;
for i=2:n-1
res_y1(i) = (y1(i+1)-2*y1(i)+y1(i-1))/dx^2 + (y1(i+1)-y1(i-1))/(2*dx) + exp(-5/Y2);
end
res_y1(n) = y1(n);
res_y2(1) = y2(1);
for i = 2:n-1
res_y2(i) = (y2(i+1)-2*y2(i)+y2(i-1))/dx^2 - y1(i).^2;
end
res_y2(n) = y2(n)-eta;
res_y3(1) = y3(1);
for i=2:n-1
res_y3(i) = (y3(i+1)-2*y3(i)+y3(i-1))/dx^2 - Y2;
end
res_y3(n) = y3(n)-1.0;
res = [res_y1;res_y2;res_y3];
end
It seems that using
res_y1(i) = (y1(i+1)-2*y1(i)+y1(i-1))/dx^2 + (y1(i+1)-y1(i-1))/(2*dx) + exp(-5/Y2);
and
res_y1(i) = (y1(i+1)-2*y1(i)+y1(i-1))/dx^2 + (y1(i+1)-y1(i-1))/(2*dx) + exp(-5/y2(i));
must be eqiuvalent, but NO. The results differ

采纳的回答

Gleb
Gleb 2021-1-26
Ok, thank you. I rewrote the code. It is more readable now. But "Unable to perform assignment because the left and right sides have a different number of elements.
Error in Gas_system_4v5/fun (line 219)
res_T_g_(ind) = (Pe_gm./Pe_g(ind)).*(T_g_(ind+1)-2*T_g_(ind)+T_g_(ind-1))/dx^2 - ...
"
again. The answer i assume is simple: cc, cg, Pe_g are matrixes n x n, not vectors, its because elementwise to n vectors gives matrix.
function Gas_system_4v5
clc
close all
n =50;
global T_fout Y_gHMXout Pe_g Pe_dk w_gHMX Pe_gm;
x0 = [800*ones(n,1);1*ones(n,1);1*ones(n,1)];
x = ((1:n)-1)/(n-1);
Pe_gm=4.8157e+004;
AlgoSelector=1;
switch AlgoSelector
case 1
Algochise='levenberg-marquardt'; % 'trust-region-dogleg'
case 2
Algochise='trust-region-dogleg' ;
case 3
Algochise='trust-region' ;
end
% options=optimoptions('fsolve', 'Display','iter', 'MaxFunctionEvaluations', 5*10^6, 'MaxIterations', 10^6, 'Algorithm', Algochise, 'StepTolerance', 10^-7, 'FunctionTolerance', 1*10^-5, 'FiniteDifferenceStepSize', 10^-6);
% options=optimoptions('fsolve', 'Display','iter', 'MaxFunctionEvaluations', 5*10^6, 'MaxIterations', 10^6, 'Algorithm', Algochise, 'StepTolerance', 10^-7, 'FunctionTolerance', 1*10^-3, 'FiniteDifferenceStepSize', 10^-4);
% 'OptimalityTolerance', 10^-11,
%options=optimoptions('fsolve', 'Display','iter', 'MaxFunctionEvaluations', 10^6, 'MaxIterations', 10^6, 'Algorithm', Algochise, 'FunctionTolerance', 10^-4, 'OptimalityTolerance', 10^-3, 'StepTolerance', 10^-5 ); % optimoptions
options=optimset( 'Display','iter', 'MaxFunEvals', 10^8, 'MaxIter', 10^6, 'Algorithm', Algochise, 'TolFun', 10^-7);
[sol,fval] = fsolve(@(x)fun(x,n),x0, options);
disp(fval);
norm(fun(sol,n))
figure
plot(x*( Pe_gm)^-0.5,sol(1:n))
hold on
figure
plot(x,sol(n+1:2*n))
hold on
figure
plot(x,sol(2*n+1:3*n))
hold on
disp(Pe_g);
disp(Pe_dk);
function res = fun(Yg,n)
x = ((1:n)-1)/(n-1);
dx = 1/(n-1);
T_fout =750;
Y_gHMXout =0.8;
cfm=5*10^4;
T_g_ = Yg(1:n);
Y_gHMX_ = Yg(n+1:2*n);
Y_HMXprod_ = Yg(2*n+1:3*n);
% for i=1:length(T_g_)
% T_g= T_g_(i);
% end
% %for i=1:n
% for i=1:length(Y_gHMX_)
% Y_gHMX= Y_gHMX_(i);
% end
%global T_fout L_g Y_gHMXout Y_B0 x Q_react_g
% format long e
Y_HMXprod0=0.2;
f1Y_cTfout=0.0;
Y_B0 = 0.0;
f_gout=0.99;
l=0.5;
f1Y_B0=Y_B0 * (1-f_gout);
selectHMX=1;
Choice1=1;
solverselect=1;
Selector=2 ;
MAX = 100;
bcinterval=[-10^2 10^2];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Y_HMXprod=0.3;
Y_gTf=0;
Y_BF3 =0;
Y_CF2=0;
Y_C=0;
Y_C2=0;
Y_cTf=0.01;
Y_B=0.01;
fg=0.7;
Q_gHMX=1359.8;
%
P =20;
m=2;
M_Tf = 100.014;
R=8.314;
ro_cTf=2.2;
M_N2O=44.013;
M_NO2=46.005;
M_NO=30.006;
M_CH2O=30.0258;
M_H2O=18.0148;
M_CO=28.01;
TmTf=603.15;
kTf=8*10^14; % C. D. Doyle, J. Appl. Polymer Sci., 5,285 (1961).
ETf=280328;
ro_B =2.34;
A=1;
%
L_g=0.1;
Y_NO2_f =0.1;
Y_NO_f=0.2;
Y_N2O_f=0.3;
Y_H2O_f=0.1;
Y_CO_f=0.1;
Y_CH2O_f=0.1;
M_B= 10.8;
M_HMX=296.1552;
M_C2F4=100;
M_BF3= 67.8;
dH_BF3 = -1135.62;
M_CF2=50;
dH_CF2=- 180; %кДж/моль
M_C=12;
dH_C=716.68 ;
M_C2=24;
dH_C2=837.74;
D = 0.5; % см2/с
wm=m/l;
dHm=cfm;
T_g=T_g_;
Y_HMXprod=Y_HMXprod_;
Y_gHMX=Y_gHMX_;
tau =@(T_g)(T_g/1000);
%_
Pn2=@(k0,k1, T_g)(k0+k1*tau(T_g));
Pn5=@(k0,k1,k2,k3,k4,T_g)(k0+k1*tau(T_g)+k2*tau(T_g).^2+k3*tau(T_g).^3+k4./tau(T_g).^2);
Ar=@(k,E)(k*exp(-E./(R*T_g)));
M_HMXprod=M_NO2*Y_NO2_f+M_NO*Y_NO_f+M_N2O*Y_N2O_f+M_H2O*Y_H2O_f+M_CO*Y_CO_f+M_CH2O*Y_CH2O_f;
Mg=(M_HMXprod*Y_HMXprod+M_HMX*Y_gHMX+M_Tf*Y_gTf + M_BF3*Y_BF3 + M_CF2*Y_CF2 +M_C*Y_C);
% плотности________________________________________________________________
rog=Mg.*P./(R.*T_g);
roc = ro_cTf*Y_cTf + ro_B*Y_B ;
ro_gf=(1-fg)*roc+fg*rog;
% теплоемкости газов______________________________________________________
c_NO2(T_g_<1200)=Pn5(16.108,75.89,-54.3874,14.30777,0.239423, T_g_(T_g_<1200))/M_NO2;
c_NO2(T_g_>=1200)=Pn5(56.82541,0.738053,-0.144721, 0.009777,-5.459911,T_g_(T_g_>=1200))/M_NO2;
c_NO(T_g_<1200)=Pn5(23.83, 12.589, -1.139, -1.497459, 0.214194, T_g_(T_g_<1200))/M_NO ;
c_NO(T_g_>=1200)=Pn5(35.9916912,0.957170,-0.148032, 0.009974,-3.004088, T_g_(T_g_>=1200))/M_NO ;
c_N2O(T_g_<1400)=(Pn5(27.67,51.148,-30.64, 6.847911, -0.157906, T_g_(T_g_<1400))/M_N2O );
c_N2O(T_g_>=1400)=Pn5(60.30274,1.034566,-0.192997, 0.012540, -6.860254, T_g_(T_g_>=1400))/M_N2O ;
c_H2O(T_g_<1700)=Pn5(30.092,6.832514,6.793435,-2.534480,0.082139,T_g_(T_g_<1700) )/M_H2O;
c_H2O(T_g_>=1700)=Pn5(41.96426,8.622053,-1.499780,0.098119,-11.15764, T_g_(T_g_>=1700) )/M_H2O;
c_CH2O(T_g_<1200)=Pn5(5.19,93.2,-44.85, 7.882279,0.551175, T_g_(T_g_<1200))/M_CH2O ;
c_CH2O(T_g_>=1200)=Pn5(71.35268,6.174497,-1.191090, 0.079564,-15.58917,T_g_(T_g_>=1200) )/M_CH2O ;
c_CO(T_g_<1300)=Pn5(25.56759,6.096130,4.054656, -2.671301, 0.131021, T_g_(T_g_<1300))/M_CO ;
c_CO(T_g_>=1300)=Pn5(35.15070,1.300095,-0.205921, 0.013550, -3.282780, T_g_(T_g_>=1300) )/M_CO ;
c_gHMX= Pn2(0.669,77.82, T_g_)/M_HMX;
C_BF3(T_g_<1000) = Pn5(21.28631, 130.3006, -109.9919, 34.28838, -0.073386, T_g_(T_g_<1000))/ M_BF3;
C_BF3(T_g_>=1000) = Pn5(81.23696, 1.096330, -0.226830, 0.015981, -6.366625, T_g_(T_g_>=1000))/M_BF3;
c_gTf(T_g_<1100)=Pn5(43.55126,175.9079,-138.7331,40.35619,-0.381260, T_g_(T_g_<1100))/M_Tf;
c_gTf(T_g_>=1100)=Pn5(129.9776,1.690918,-0.340087,0.023448,-10.83204, T_g_(T_g_>=1100))/M_Tf;
C_CF2(T_g_<600) =Pn5(8.825772, 125.3652, -129.4940, 50.33101, 0.259749, T_g_(T_g_<600))/ M_CF2;
C_CF2(T_g_>=600) = Pn5(59.34753, -2.317176, 0.890518, -0.055879, -3.467545, T_g_(T_g_>=600))/M_CF2;
C_C= Pn5(21.17510, -0.812428, 0.448537, -0.043256, -0.013103, T_g_)/M_C;
C_C2(T_g_<700) = Pn5(123.9092,-348.0067,485.0971, -232.7994, -1.240298, T_g_(T_g_<700))/M_C2;
C_C2(T_g_>=700) = Pn5(30.50408,5.445811,-0.853373, 0.065641, 0.814750, T_g_(T_g_>=700))/M_C2;
C_HMXproducts=c_NO2*Y_NO2_f+c_NO*Y_NO_f+c_N2O*Y_N2O_f+c_H2O*Y_H2O_f+c_CO*Y_CO_f+c_CH2O*Y_CH2O_f;
cg=C_HMXproducts.*Y_HMXprod+c_gHMX.*Y_gHMX+c_gTf.*Y_gTf+C_BF3.*Y_BF3 + C_CF2.*Y_CF2 + C_C.*Y_C+C_C2.*Y_C2;
% теплоемкости к-веществ__________________________________________________
c_cTf(T_g_<TmTf)=1.04;
c_cTf(T_g_>=TmTf)= (0.61488+0.001194*1000*tau(T_g));
c_B = (10.18574 + 29.24415*tau(T_g) -18.02137*tau(T_g).^2 +4.212326*tau(T_g).^3 )/10.8 ;
cc=c_cTf.*Y_cTf+ c_B.*Y_B;
c_gf=((1-fg).*roc.*cc+fg.*rog.*cg)./ro_gf;
% теплопроводности_________________________________________________________
lc=2.5e-3*Y_cTf + 0.27*Y_B;
l_HMXproducts = 2e-4*Y_NO2_f+2.59e-4*Y_NO_f+1.74e-4*Y_N2O_f+2.79e-4*Y_H2O_f+2.5e-4*Y_CO_f+ 2e-4*Y_CH2O_f; % l_BF3
lg= l_HMXproducts*Y_HMXprod+12.5e-4*Y_gHMX+2.5e-4*Y_gTf +2e-4*(Y_BF3 +Y_C+Y_C2 +Y_CF2);
l_gf= ((1-fg).*lc.*cc+fg.*lg.*cg)./((1-fg).*cc+fg.*cg);
%________________________________________________________________________
% кинетмка
%Tf->C2F4
% C2F4->2CF2
% B+3CF2->2BF3+3C
% 2C->C2
switch selectHMX
case 1
G_gHMX = fg.*ro_gf.*Y_gHMX.*Ar(10^14.2, 165268)/M_HMX; %
case 2
G_gHMX =fg.*ro_gf.*Y_gHMX.*Ar(10^12.5, 158992)/M_HMX;
case 3
G_gHMX =fg.*ro_gf.*Y_gHMX.*Ar(5*10^13, 158992)/M_HMX;
end
G_Tf= ((1-fg).*ro_cTf.*Y_cTf.*Ar(kTf,ETf)/M_Tf);
G_CF2dec = fg.*ro_gf.*Y_gTf.*T_g.^(-0.5).*Ar(7.82*10^15,2.33*10^5)/M_C2F4;
G_B_CF2 = fg*ro_gf.*Y_CF2.*(Y_B*ro_B).^3.*T_g.^(0.5).*Ar(4*10^14,8.37*10^4)/(M_CF2*M_B^3);
G_C = fg.*ro_gf.^2.*Y_C.^2.*T_g.^(-1.6).*Ar(1.80*10^21,0)/M_C;
w_gHMX= -M_HMX*G_gHMX;
w_gHMXprod = M_HMX*G_gHMX;
%Tf->C2F4
w_cTf=-M_Tf*G_Tf;
w_gTf =-M_Tf*w_cTf;
% B+3CF2->BF3+3C
w_B=-M_B*G_B_CF2;
% C2F4->2CF2, B+3CF2->2BF3+3C
w_CF2=2*M_CF2*G_CF2dec-3*M_CF2*G_B_CF2;
% B+3CF2->2BF3+3C
w_BF3 =2*M_BF3* G_B_CF2;
% B+3CF2->2BF3+3C, 2C->C2
w_C = -2*M_C*G_C +3*M_C*G_B_CF2;
w_C2 = M_C2*G_C;
G_condphase= -w_cTf-w_B;
% тепловыделение___________________________________________________________
Q_react_g= Q_gHMX*w_gHMX/(wm*dHm); %0*( -w_cTf*824/M_Tf - w_gTf*658.56/M_Tf + w_CF2*dH_CF2/M_CF2 + w_BF3*dH_BF3/M_BF3 + dH_C*w_C/M_C +dH_C2*w_C2/M_C2)*1e3/(wm*dHm)));
%disp(Q_react_g);
Pe_g= cfm*m.*l./l_gf;
Pe_dk = m*l./(D*rog);
%________________________________________________________________________
res_T_g_ = zeros(n,1);
res_Y_gHMX_ = zeros(n,1);
res_Y_HMXprod_=zeros(n,1);
%________________________________________________________________________
ind = 2:n-1;
res_T_g_(1) = T_g_(1)-T_fout;
res_T_g_(ind) = (Pe_gm./Pe_g(ind)).*(T_g_(ind+1)-2*T_g_(ind)+T_g_(ind-1))/dx^2 - ...
(Pe_gm)^0.5 .* (c_gf(ind)/cfm).*(T_g_(ind+1)-T_g_(ind))/(2*dx) - Q_react_g(ind);
res_T_g_(n) = ( T_g_(n)-T_g_(n-1))/dx ;
%________________________________________________________________________
ind = 2:n-1;
res_Y_gHMX_(1) = Y_gHMX_(1)-Y_gHMXout;
res_Y_gHMX_(ind) =(1./Pe_dk(ind))*(Y_gHMX_(ind+1)-2*Y_gHMX_(ind)+Y_gHMX_(ind-1))/dx^2-(Y_gHMX_(ind+1)-Y_gHMX_(ind-1))/(2*dx)+w_gHMX(ind)/wm; %w_gHMX Y_gHMX_(i).*ro_gf.*(10^14.2)*exp(-165268./(8.314*T_g_(i) ))
res_Y_gHMX_(n) =( Y_gHMX_(n)-Y_gHMX_(n-1))/dx;
%________________________________________________________________________
ind = 2:n-1;
res_Y_HMXprod_(1) = Y_HMXprod_(1)-Y_HMXprod0;
res_Y_HMXprod_(ind) =(1./Pe_dk(ind))*(Y_HMXprod_(ind+1)-2*Y_HMXprod_(ind)+Y_HMXprod_(ind-1))/dx^2-(Y_HMXprod_(ind+1)-Y_HMXprod_(ind-1))/(2*dx)+w_gHMXprod(ind)/wm; %w_gHMX Y_gHMX_(i).*ro_gf.*(10^14.2)*exp(-165268./(8.314*T_g_(i) ))
res_Y_HMXprod_(n) =( Y_HMXprod_(n)-Y_HMXprod_(n-1))/dx;
res = [res_T_g_;res_Y_gHMX_; res_Y_HMXprod_;];
end
end
  4 个评论
Cris LaPierre
Cris LaPierre 2021-1-28
I put this in the answer: "I also had to be sure all vectors were nx1 to avoid this issue, so I assign using c_NO2(T_g_<1200,1)=..."
This syntax places the returned values into a single column instead of a row. This was necessary to prevent implicit explansion turning the results of some calculations into 50x50 matrices.
As for the second question, this was how I elected to make the 2 variables returned by the function (Pe_g, Pe_dk) available for the disp commands (I commented them out). You used a global variable, but it's best to avoid globals if at all possible.
If you remove the disp commands, you can remove [~,Pe_g, Pe_dk] = fun(sol,n) as well.
Gleb
Gleb 2021-2-7
Thank you. Problem solved. Regards.

请先登录,再进行评论。

更多回答(1 个)

Cris LaPierre
Cris LaPierre 2021-1-20
Why would you expect them to be the same?
When you get around to calculation res_y1(i), Y2 has the value of y2(length(y1)). If you substitute that in, you'll see the last part fo the equation is very different
exp(-5/y2(length(y1)));
%vs
exp(-5/y2(i));
The first uses the same value of Y2 for every loop, while the second uses a different element of y2 for each loop.
  26 个评论
Cris LaPierre
Cris LaPierre 2021-1-26
编辑:Cris LaPierre 2021-1-26
I've shown you how to remove all the for loops, so I think you only need 4 anonymous functions in your code.
tau =@(T_g)(T_g/1000);
%_
Pn2=@(k0,k1, T_g)(k0+k1*tau(T_g));
Pn5=@(k0,k1,k2,k3,k4, T_g)(k0+k1*tau(T_g)+k2*tau(T_g).^2+k3*tau(T_g).^3+k4./tau(T_g).^2);
Ar=@(k,E)(k*exp(-E./(R*T_g_)));
See if you can rewrite your code to remove all other anonymous functions.
Cris LaPierre
Cris LaPierre 2021-1-26
is there some advantage over traditional for loop?
Loops = slower code. In MATLAB, you can often avoid them by taking advantage of elementwise operators.
Well... There are 3 equations, and 3 loops, isnt it?
I've got to leave something for you to do. You can use the example I shared to see how to remove the loops from the remaining 2 equations.

请先登录,再进行评论。

标签

Community Treasure Hunt

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

Start Hunting!

Translated by