problem in the for - loop

2 Ansichten (letzte 30 Tage)
Gleb
Gleb am 20 Jan. 2021
Kommentiert: Gleb am 7 Feb. 2021
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

Akzeptierte Antwort

Gleb
Gleb am 26 Jan. 2021
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 Kommentare
Cris LaPierre
Cris LaPierre am 28 Jan. 2021
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 am 7 Feb. 2021
Thank you. Problem solved. Regards.

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (1)

Cris LaPierre
Cris LaPierre am 20 Jan. 2021
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 Kommentare
Cris LaPierre
Cris LaPierre am 26 Jan. 2021
Bearbeitet: Cris LaPierre am 26 Jan. 2021
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 am 26 Jan. 2021
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.

Melden Sie sich an, um zu kommentieren.

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by