What can I change on the code below to make it run faster?

조회 수: 1 (최근 30일)
Dursman Mchabe
Dursman Mchabe 2019년 2월 1일
댓글: Steven Lord 2019년 2월 1일
Hello everyone,
What can I change on the code below to make it run faster?
function SlurryCaseODE45Jan31
function C=kinetics(theta,t)
c0=[0.076298787;1e-23;1e-23;1e-23;1e-23;47.61362376;1e-23];
[T,Cv]=ode45(@DifEq,t,c0);
%
function dC=DifEq(t,c)
dcdt=zeros(7,1);
dcdt(1) = (1 /V_Headspace).* (F.* CSO2_in - F.* c(1) ) - ((c(1).* R.* Temp - HSO2.* ((c(3).* c(7).^2)/(c(7).^2 + KSO2.* c(7) + KSO2.* KHSO3))) / ((1/theta(1)) + HSO2/ ((1 + DCa2.* c(5) / (DSO2.* c(3))).* theta(2)) ));
dcdt(2) = (1 /V_Headspace).* (F .* CCO2_in - F .* c(2)) - ((theta(3).* (1 + DCa2 * c(5) / (DCO2 .* c(4)))/ HCO2).* (c(2) .* R .* Temp - HCO2.*((c(4).* c(7).^2)/(c(7).^2 + KCO2.* c(7) + KCO2.* KHCO3))));
dcdt(3) = ((c(1).* R.* Temp - HSO2.* ((c(3).* c(7).^2)/(c(7).^2 + KSO2.* c(7) + KSO2.* KHSO3))) / ((1/theta(1)) + HSO2/ ((1 + DCa2.* c(5) / (DSO2.* c(3))).* theta(2)) )) - (theta(10) .*exp(-theta(11)/Temp).*((c(5).*((c(3).* KSO2.* KHSO3)/(c(7).^2 + KSO2 .* c(7) + KSO2 .*KHSO3))/theta(4)) - 1).^theta(9) .*(( theta(5) .* theta(4))/(c(5).*((c(3).* KSO2.* KHSO3)/(c(7).^2 + KSO2 .* c(7) + KSO2 .*KHSO3)))));
dcdt(4) = (-theta(6).*theta(7).*MWCaCO3.*c(6).*(- 2.*c(5) + ((c(3).* KSO2.* c(7))/(c(7).^2 + KSO2 .* c(7) + KSO2 .*KHSO3)) + 2.*((c(3).* KSO2.* KHSO3)/(c(7).^2 + KSO2 .* c(7) + KSO2 .*KHSO3)) + ((c(4).* KCO2.* c(7))/(c(7).^2 + KCO2 .* c(7) + KCO2 .*KHCO3)) + 2.*((c(4).* KCO2.* KHCO3)/(c(7).^2 + KCO2 .* c(7) + KCO2 .*KHCO3)) + (Kw/c(7))).*(1-(theta(8).*(- 2.*c(5) + ((c(3).* KSO2.* c(7))/(c(7).^2 + KSO2 .* c(7) + KSO2 .*KHSO3)) + 2.*((c(3).* KSO2.* KHSO3)/(c(7).^2 + KSO2 .* c(7) + KSO2 .*KHSO3)) + ((c(4).* KCO2.* c(7))/(c(7).^2 + KCO2 .* c(7) + KCO2 .*KHCO3)) + 2.*((c(4).* KCO2.* KHCO3)/(c(7).^2 + KCO2 .* c(7) + KCO2 .*KHCO3)) + (Kw/c(7))))/(1+ theta(8).*(- 2.*c(5) + ((c(3).* KSO2.* c(7))/(c(7).^2 + KSO2 .* c(7) + KSO2 .*KHSO3)) + 2.*((c(3).* KSO2.* KHSO3)/(c(7).^2 + KSO2 .* c(7) + KSO2 .*KHSO3)) + ((c(4).* KCO2.* c(7))/(c(7).^2 + KCO2 .* c(7) + KCO2 .*KHCO3)) + 2.*((c(4).* KCO2.* KHCO3)/(c(7).^2 + KCO2 .* c(7) + KCO2 .*KHCO3)) + (Kw/c(7)))))) - ((theta(3).* (1 + DCa2 * c(5) / (DCO2 .* c(4)))/ HCO2).* (c(2) .* R .* Temp - HCO2.*((c(4).* c(7).^2)/(c(7).^2 + KCO2.* c(7) + KCO2.* KHCO3))));
dcdt(5) = (-theta(6).*theta(7).*MWCaCO3.*c(6).*(- 2.*c(5) + ((c(3).* KSO2.* c(7))/(c(7).^2 + KSO2 .* c(7) + KSO2 .*KHSO3)) + 2.*((c(3).* KSO2.* KHSO3)/(c(7).^2 + KSO2 .* c(7) + KSO2 .*KHSO3)) + ((c(4).* KCO2.* c(7))/(c(7).^2 + KCO2 .* c(7) + KCO2 .*KHCO3)) + 2.*((c(4).* KCO2.* KHCO3)/(c(7).^2 + KCO2 .* c(7) + KCO2 .*KHCO3)) + (Kw/c(7))).*(1-(theta(8).*(- 2.*c(5) + ((c(3).* KSO2.* c(7))/(c(7).^2 + KSO2 .* c(7) + KSO2 .*KHSO3)) + 2.*((c(3).* KSO2.* KHSO3)/(c(7).^2 + KSO2 .* c(7) + KSO2 .*KHSO3)) + ((c(4).* KCO2.* c(7))/(c(7).^2 + KCO2 .* c(7) + KCO2 .*KHCO3)) + 2.*((c(4).* KCO2.* KHCO3)/(c(7).^2 + KCO2 .* c(7) + KCO2 .*KHCO3)) + (Kw/c(7))))/(1+ theta(8).*(- 2.*c(5) + ((c(3).* KSO2.* c(7))/(c(7).^2 + KSO2 .* c(7) + KSO2 .*KHSO3)) + 2.*((c(3).* KSO2.* KHSO3)/(c(7).^2 + KSO2 .* c(7) + KSO2 .*KHSO3)) + ((c(4).* KCO2.* c(7))/(c(7).^2 + KCO2 .* c(7) + KCO2 .*KHCO3)) + 2.*((c(4).* KCO2.* KHCO3)/(c(7).^2 + KCO2 .* c(7) + KCO2 .*KHCO3)) + (Kw/c(7)))))) - (theta(10) .*exp(-theta(11)/Temp).*((c(5).*((c(3).* KSO2.* KHSO3)/(c(7).^2 + KSO2 .* c(7) + KSO2 .*KHSO3))/theta(4)) - 1).^theta(9) .*(( theta(5) .* theta(4))/(c(5).*((c(3).* KSO2.* KHSO3)/(c(7).^2 + KSO2 .* c(7) + KSO2 .*KHSO3)))));
dcdt(6) = (-theta(6).*theta(7).*MWCaCO3.*c(6).*(- 2.*c(5) + ((c(3).* KSO2.* c(7))/(c(7).^2 + KSO2 .* c(7) + KSO2 .*KHSO3)) + 2.*((c(3).* KSO2.* KHSO3)/(c(7).^2 + KSO2 .* c(7) + KSO2 .*KHSO3)) + ((c(4).* KCO2.* c(7))/(c(7).^2 + KCO2 .* c(7) + KCO2 .*KHCO3)) + 2.*((c(4).* KCO2.* KHCO3)/(c(7).^2 + KCO2 .* c(7) + KCO2 .*KHCO3)) + (Kw/c(7))).*(1-(theta(8).*(- 2.*c(5) + ((c(3).* KSO2.* c(7))/(c(7).^2 + KSO2 .* c(7) + KSO2 .*KHSO3)) + 2.*((c(3).* KSO2.* KHSO3)/(c(7).^2 + KSO2 .* c(7) + KSO2 .*KHSO3)) + ((c(4).* KCO2.* c(7))/(c(7).^2 + KCO2 .* c(7) + KCO2 .*KHCO3)) + 2.*((c(4).* KCO2.* KHCO3)/(c(7).^2 + KCO2 .* c(7) + KCO2 .*KHCO3)) + (Kw/c(7))))/(1+ theta(8).*(- 2.*c(5) + ((c(3).* KSO2.* c(7))/(c(7).^2 + KSO2 .* c(7) + KSO2 .*KHSO3)) + 2.*((c(3).* KSO2.* KHSO3)/(c(7).^2 + KSO2 .* c(7) + KSO2 .*KHSO3)) + ((c(4).* KCO2.* c(7))/(c(7).^2 + KCO2 .* c(7) + KCO2 .*KHCO3)) + 2.*((c(4).* KCO2.* KHCO3)/(c(7).^2 + KCO2 .* c(7) + KCO2 .*KHCO3)) + (Kw/c(7))))));
dcdt(7) = (theta(10) .*exp(-theta(11)/Temp).*((c(5).*((c(3).* KSO2.* KHSO3)/(c(7).^2 + KSO2 .* c(7) + KSO2 .*KHSO3))/theta(4)) - 1).^theta(9) .*(( theta(5) .* theta(4))/(c(5).*((c(3).* KSO2.* KHSO3)/(c(7).^2 + KSO2 .* c(7) + KSO2 .*KHSO3)))));
dC=dcdt;
end
C=Cv;
end
t=[0
6000
12000
18000
24000
30000
36000
42000
48000
54000
60000
66000
72000
78000
84000
90000
96000
102000
108000
114000
120000
126000
132000
138000
144000
150000
156000
162000
168000
174000
180000];
c=[0.076298787 0 0 0 0 47.61362376 0
0.019461171 0.029472654 6.668799696 16.23121199 9.959999844 35.99050802 6.300684797
0.020167931 0.042887117 11.84949115 24.04825592 14.75631477 29.46171885 9.33482827
0.019107791 0.051014128 17.25229685 29.01821217 17.80549194 25.47239217 11.26373435
0.019107791 0.051014128 17.25229685 29.01821217 17.80549194 25.47239217 11.26373435
0.019343503 0.063743181 22.80780219 36.81976241 22.59189684 18.84582265 14.29160875
0.019343503 0.078528466 28.5615147 46.11097327 28.29200863 10.86857406 17.89749311
0.019696506 0.081955518 34.59038619 49.09849644 30.12402642 9.128367823 19.05642553
0.020850931 0.088613792 40.44036333 53.76834115 32.98847955 5.388343709 20.86847538
0.024902168 0.092530424 46.03725425 56.87434597 34.89336965 3.160034015 22.07350674
0.024902168 0.092530424 46.03725425 56.87434597 34.89336965 3.160034015 22.07350674
0.031049924 0.097915792 51.29729887 61.13770121 37.50806004 0 23.72755696
0.035572586 0.097915792 56.06755651 61.90018578 37.97509734 0 24.02300424
0.038601396 0.097915792 60.70353205 62.73848974 38.48857555 0 24.34782998
0.042662438 0.097915792 65.10035551 63.65803386 39.05181501 0 24.70413464
0.047985767 0.097915792 68.80560315 64.36554967 39.48518283 0 24.97828263
0.047985767 0.097915792 68.80560315 64.36554967 39.48518283 0 24.97828263
0.051943095 0.097915792 72.44446145 65.33369776 40.07819321 0 25.35342034
0.056347712 0.097915792 76.00361827 66.58555956 40.84498401 0 25.83849134
0.060163613 0.097915792 78.89271653 67.55634934 41.43961248 0 26.21465265
0.063414181 0.097915792 81.66324156 68.6913213 42.13480587 0 26.65443122
0.06666475 0.097915792 84.24007982 69.93531357 42.89677643 0 27.13645294
0.06666475 0.097915792 84.24007982 69.93531357 42.89677643 0 27.13645294
0.069208558 0.097915792 86.52721476 71.15222604 43.64216005 0 27.60798179
0.071540791 0.097915792 88.61843725 72.41217703 44.41390563 0 28.09618718
0.073377839 0.097915792 90.46782755 73.63940812 45.16560962 0 28.57171429
0.074673315 0.097915792 92.12354968 74.82827386 45.89381405 0 29.0323756
0.07512098 0.097915792 92.26733175 74.82827386 45.89381405 0 29.0323756
0.07512098 0.097915792 92.26733175 74.82827386 45.89381405 0 29.0323756
0.075709695 0.097915792 92.33853494 74.82827386 45.89381405 0 29.0323756
0.076063075 0.097915792 92.36679575 74.82827386 45.89381405 0 29.0323756];
% theta0 = [ 1.95e-3; 1.706e-5; 7.28e-2; 1e-23; 8.314; 323.15; 149; 4.23e-6; 1.39e-9; 2.89e-9; 5.78e-3; 9.598e-4; 5.15e3; 3.53e-9; 1.07e-7; 10; 6.387e-3; 12.54; 100.0869; 1.03493; 1.7e-3; 6.55e-8; 6.24; 5.68e-5; 5.3e-8;3];
theta0 = [4.23e-6;5.78e-3;9.598e-4;1.07e-7;10;6.387e-3;12.54;1.03493;3;0.162;5153];
V_Headspace = 1.95e-3;
F = 1.706e-5;
CSO2_in = 7.28e-2;
CCO2_in = 1e-23;
R = 8.314;
Temp = 323.15;
HSO2 = 149;
%theta(1) = 4.23e-6; % kga
DCa2 = 1.39e-9;
DSO2 = 2.89e-9;
%theta(2) = 5.78e-3; % kLa_SO2
%theta(3) = 9.598e-4; % kLa_CO2
HCO2 = 5.15e3;
DCO2 = 3.53e-9;
%theta(4) = 1.07e-7; %KspCaSO3
%theta(5) = 10; % BETCaSO3
%theta(6) = 6.387e-3; % ktot
%theta(7) = 12.54; % BETCaCO3
MWCaCO3 = 100.0869;
%theta(8) = 1.03493; % Kad
KCO2 = 1.7e-3;
KHCO3 = 6.55e-8;
KSO2 = 6.24;
KHSO3 = 5.68e-5;
Kw = 5.3e-8;
%theta(9) = 3; % n
%theta(10) = 0.162; % A0
%theta(11) = 5153; % Ea
%theta0 = [1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1];
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);
fprintf(1,'\tConstants:\n')
for k1 = 1:length(theta)
fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
end
tv = linspace(min(t), max(t),31);
%tv = t;
Cfit = kinetics(theta, tv);
figure(1)
plot(t, c, 'p')
hold on
hlp = plot(tv, Cfit);
hold off
grid
xlabel('Time')
ylabel('Concentration')
end

이 질문에 답변하려면 로그인하십시오.

채택된 답변

Jan
Jan 2019년 2월 1일
편집: Jan 2019년 2월 1일
A small improvement is possible, when you exploit, that the terms of dcdt(4:6) are identical, but I assume this saves some percent only. A better improvement can be achieved, if you define the tolerances of ODE45 exactly to your needs instead of using the default values. This influences the number of integration steps remarkably and can reduce the runtime substantially.
But before you start to optimize the ODE45 call, use the profiler to find out, if this is the bottleneck at all: Care only for the part, which consumes the most time. If a piece of code needs 2% of the total time only, reducing its runtime by 50% will let the program run 1% faster only.
So please mention, what profiler reveals.
  댓글 수: 2
Dursman Mchabe
Dursman Mchabe 2019년 2월 1일
Thanks a lot Jan. I trust that this suggestions will help.
Steven Lord
Steven Lord 2019년 2월 1일
Another possibility would be to try a stiffer solver if the ODE solving is the bottleneck. I'm not sure how stiff your equations are, but it's an easy technique to try. Just change the function you're calling from ode45 to ode15s.
If the ODE solver itself isn't the bottleneck but lsqcurvefit calling the ODE solver over and over is, you might want to adjust some of the lsqcurvefit options.

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

카테고리

Help Center 및 File Exchange에서 Performance and Memory에 대해 자세히 알아보기

태그

Community Treasure Hunt

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

Start Hunting!

Translated by