problem in the for - loop
Show older comments
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
Accepted Answer
Gleb
on 26 Jan 2021
4 Comments
You are correct. i had to make sure to use the elementwise operator (.* and ./) in calculating those variables. I also had to be sure all vectors were nx1 to avoid this issue, so I assign using c_NO2(T_g_<1200,1)=...
Here's what my solution looks like. Notice I also avoid using global variables. I unfortunately have no idea if the result is correct.
Gas_system_4v3
function Gas_system_4v3
n =50;
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);
% Call the function with the solution values to capture Pe_g and Pe_dk
[~,Pe_g, Pe_dk] = fun(sol,n);
% 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, Pe_g, Pe_dk] = 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
% 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;
Y_B=0;
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;
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,1)=Pn5(16.108,75.89,-54.3874,14.30777,0.239423,T_g_(T_g_<1200))/M_NO2;
c_NO2(T_g_>=1200,1)=Pn5(56.82541,0.738053,-0.144721, 0.009777,-5.459911,T_g_(T_g_>=1200))/M_NO2;
c_NO(T_g_<1200,1)=Pn5(23.83, 12.589, -1.139, -1.497459, 0.214194,T_g_(T_g_<1200))/M_NO;
c_NO(T_g_>=1200,1)=Pn5(35.9916912,0.957170,-0.148032, 0.009974,-3.004088,T_g_(T_g_>=1200))/M_NO;
c_N2O(T_g_<1400,1)=Pn5(27.67,51.148,-30.64, 6.847911, -0.157906,T_g_(T_g_<1400))/M_N2O;
c_N2O(T_g_>=1400,1)=Pn5(60.30274,1.034566,-0.192997, 0.012540, -6.860254,T_g_(T_g_>=1400))/M_N2O;
c_H2O(T_g_<1700,1)=Pn5(30.092,6.832514,6.793435,-2.534480,0.082139,T_g_(T_g_<1700))/M_H2O;
c_H2O(T_g_>=1700,1)=Pn5(41.96426,8.622053,-1.499780,0.098119,-11.15764,T_g_(T_g_>=1700))/M_H2O;
c_CH2O(T_g_<1200,1)=Pn5(5.19,93.2,-44.85, 7.882279,0.551175,T_g_(T_g_<1200))/M_CH2O;
c_CH2O(T_g_>=1200,1)=Pn5(71.35268,6.174497,-1.191090, 0.079564,-15.58917,T_g_(T_g_>=1200))/M_CH2O;
c_CO(T_g_<1300,1)=Pn5(25.56759,6.096130,4.054656, -2.671301, 0.131021,T_g_(T_g_<1300))/M_CO;
c_CO(T_g_>=1300,1)=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,1) = Pn5(21.28631, 130.3006, -109.9919, 34.28838, -0.073386,T_g_(T_g_<1000))/ M_BF3;
C_BF3(T_g_>=1000,1) = Pn5(81.23696, 1.096330, -0.226830, 0.015981, -6.366625,T_g_(T_g_>=1000))/M_BF3;
c_gTf(T_g_<1100,1)=Pn5(43.55126,175.9079,-138.7331,40.35619,-0.381260,T_g_(T_g_<1100))/M_Tf;
c_gTf(T_g_>=1100,1)=Pn5(129.9776,1.690918,-0.340087,0.023448,-10.83204,T_g_(T_g_>=1100))/M_Tf;
C_CF2(T_g_<600,1) = Pn5(8.825772, 125.3652, -129.4940, 50.33101, 0.259749,T_g_(T_g_<600))/ M_CF2;
C_CF2(T_g_>=600,1) = 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,1) = Pn5(123.9092,-348.0067,485.0971, -232.7994, -1.240298,T_g_(T_g_<700))/M_C2;
C_C2(T_g_>=700,1) = 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)=1.04;
c_cTf(T_g_>=TmTf,1)= (0.61488+0.001194*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+ ... $ -258
+ 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 + 0*(1./Pe_dk(ind)).*(c_gf(ind)/cfm).*ro_gf(ind).*((T_g_(ind+1)-T_g_(ind-1))/(2*dx)).*Y_gHMX_(ind).*(Y_gHMX_(ind+1)-Y_gHMX_(ind-1))/(2*dx) - ...
(Pe_gm)^0.5 .* (c_gf(ind)/cfm).*(T_g_(ind+1)-T_g_(ind))/(2*dx) - Q_react_g(ind); % Q_gHMX*Y_gHMX_(i).*ro_gf.*Ar(10^14.2, 165268) 1*
res_T_g_(n) = ( T_g_(n)-T_g_(n-1))/dx ;
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;
res_Y_HMXprod_(1) = Y_HMXprod_(1)-Y_HMXprod0;
% w_gHMXprod= fg*Y_gHMX_(i).*ro_gf.*(10^14.2)*exp(-165268./(8.314*T_g_(i) ));
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
Gleb
on 27 Jan 2021
Cris LaPierre
on 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
on 7 Feb 2021
More Answers (1)
Cris LaPierre
on 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 Comments
Cris LaPierre
on 20 Jan 2021
The underlying problem appears to still be the same. In one case, Y2 is always the same value for every loop while in the other, y2(i) will be a different value in every loop.
Gleb
on 21 Jan 2021
Cris LaPierre
on 21 Jan 2021
You can absolutely use a function inside a for loop. You can't define an in-file function in the loop, but that's a different issue.
I guess it's all in how you define your anonymous function. If you define your anonymous function in a variable, you can use it as a function elsewhere in your code.
% define an anonymous function
sqr = @(x) x.^2;
% use the anonymous function
a = sqr(5)
Gleb
on 25 Jan 2021
Cris LaPierre
on 25 Jan 2021
No, I don't think I would code it like this. Of course, you have left out a lot of the details, so it's hard to say for certain. Generally, I wouldn't define the anonymous function inside a for loop.
Also, you are not following the syntax shared in the example. You should indicate a variable used in the equation right after the @, then use that variable in the function. When you call the function, that is when you use your actual indexed variable.
So based on your example, I think it would look like this (untested).
res_Y1(1) = Y1(1)-Y0;
A=@(myY1, myY2) (yY1+myY); % some expressions
B=@(myY1, myY2) (yY1.*myY);
Q=@(myY1, myY2) A(myY1, myY2)+B(myY1, myY2);
for i = 2:n-1
res_Y1(i) =(Y1(i+1)-2*Y1(i)+Y1(i-1))/dx^2-Q(Y1(i), Y2(i))
end
res_Y1(n) =( Y1(n)-Y1(n-1))/dx;
% loop 2
res_Y2(1) = Y2(1)-Y20;
for i = 2:n-1
res_Y2(i) =(Y2(i+1)-2*Y2(i)+Y2(i-1))/dx^2+Q(Y1(i), Y2(i))
end
...
Gleb
on 25 Jan 2021
Cris LaPierre
on 25 Jan 2021
The variable myY1 is a variable that only exists in the anonymous function. When you call A, that is when you pass in the actual value.
Also note that here you have only defined functions. Be sure to also use them to perform an actual calculation.
if Y1(i)<1000
% define your function
A=@(myY1, myY2) (myY1+myY2^2); % some expressions
% use your function
res(i) = A(Y1(i),Y2(i));
else
% define your function
A=@(myY1, myY2) (2*myY1+3*myY2^2)
% use your function
res(i) = A(Y1(i),Y2(i));
end
Gleb
on 25 Jan 2021
Cris LaPierre
on 25 Jan 2021
For that matter, Y1(i) is never defined in the code you shared. This is the problem with trying to simplify your example. If you want an answer that is guaranteed to work, you need to share your actual code.
Cris LaPierre
on 25 Jan 2021
Edited: Cris LaPierre
on 26 Jan 2021
I see a few things.
- First, your if statement uses T_g_(i), but i has not been defined.
- You are defining a new function, c_NO2. However, you try to use it as a variable.
- c_NO2 calls Pn5 with 5 inputs, but you defined it as having 6.
I suggest reviewing the link I shared about function basics. In particular, how to call a function. There are several anonymous functions created in the code that are not called properly.
- tau (no input provided in Pn5)
- Pn2 (has 3 inputs, but you call it with 2)
- Pn5 (has 6 inputs, but you call it with 5)
- Ar (has 3 inputs, but you call it with 2)
- ro_gf (has 3 inputs, but you often use it without any inputs
One more thing to figure out is what is T_g. You use it in your if statements conditional, but have never defined it. It is a variable used in some of your anonymous functions, but it only exists within the function's scope.
After fixing these issues, I think my approach to what you have in line 156 would use logical indexing instead. This allows me to compute a value for all values that meet the condition at once, rather than having to loop through each value.
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;
Notice that c_NO2 is a variable that captures the result of the calculation performed by the function Pn5.
Gleb
on 25 Jan 2021
Cris LaPierre
on 26 Jan 2021
cc is a variable that got incorrectly set up as an anonymous function. I think it should be
cc=c_cTf.*Y_cTf+ c_B(T_g_).*Y_B;
but you need to confirm that.
Gleb
on 26 Jan 2021
Cris LaPierre
on 26 Jan 2021
Edited: Cris LaPierre
on 26 Jan 2021
Not when I run it.
Check the dimensions of your variables to identify the differences. Then go back to where each of these variables were created to find any errors.
Gleb
on 26 Jan 2021
Cris LaPierre
on 26 Jan 2021
Maybe a new error, but same issue. I think c_gf should not be an anonymous function. I think you are actually trying to create a vector of values.
Cris LaPierre
on 26 Jan 2021
Edited: Cris LaPierre
on 26 Jan 2021
Perhaps a better way to approach this is start with the first anonymous functions your create (lines 135-148).
dHm=cfm;
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, T_g)(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=@(T_g, Y_HMXprod, Y_gHMX)(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=@(T_g, Y_HMXprod, Y_gHMX)(Mg(T_g, Y_HMXprod, Y_gHMX).*P./(R.*T_g));
roc = ro_cTf*Y_cTf + ro_B*Y_B ;
ro_gf=@(T_g, Y_HMXprod, Y_gHMX)((1-fg)*roc+fg*rog(T_g, Y_HMXprod, Y_gHMX));
We've already talked about tau, Pn2 and Pn5. They should remain anonymous functions. Take a hard look at the rest. Which do you think should be vectors of values instead of anonymous functions?
A good criteria to consider is if you ever create the inputs used when you call the function, or if you call the function multiple times with different input values.
Cris LaPierre
on 26 Jan 2021
Edited: Cris LaPierre
on 26 Jan 2021
Let's focus on the ones in the code I shared for now.
As for the for loops at the bottom, you don't need them. I would use indexing.
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 + 0*(1./Pe_dk(ind)).*(c_gf(ind)/cfm).*ro_gf(ind).*((T_g_(ind+1)-T_g_(ind-1))/(2*dx)).*Y_gHMX_(ind).*(Y_gHMX_(ind+1)-Y_gHMX_(ind-1))/(2*dx) - ...
(Pe_gm)^0.5 .* (c_gf(ind)/cfm).*(T_g_(ind+1)-T_g_(ind))/(2*dx) - Q_react_g(ind);
I couldn't help but notice that a large section of your calculation of res_T_g_ is multiplied by 0. Since this doesn't contribute anything to the calculation, you could simply the code by removing it.
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);
Gleb
on 26 Jan 2021
Cris LaPierre
on 26 Jan 2021
Edited: Cris LaPierre
on 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
on 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.
Categories
Find more on Linear Least Squares in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
