Iteration over initial condition in ode45
显示 更早的评论
Hi, I am trying to solve the following differential equation using ode45, which works fine but I need to iterate the initial condition of ti for different values of e (0.4,0.5,0.6,0.7,0.8,0.9) and Pg (3,4,5) untile the solution reach the value of 0.05. For instance for the case of e = 0.73 and Pg = 3, only with ti = 984, the solution will reach 0.05. My question is how can I solve this eqaution iteratively until for each values of e and Pg it finds the ti that gives the CDF = 0.05.
I know I asked this question before here but still couldn't figure out the solution for it. I would appreciate if anyone can help me. Thanks!
clc; clear;
ti = 984;
CDF_in = 0;
tEnd = 7200;
[tsol, CDFsol] = ode45(@(t,CDF) firstODEfun(t,CDF), [ti tEnd], CDF_in);
plot(tsol, CDFsol)
function dCDF = firstODEfun(t,~)
t0 = 2.512E-018;
r0 = 8.275E-02;
Nrod = 91;
l = 1.59E-02;
drod = 0.0135;
Drod = 0.01174;
sigma = 5.67E-08;
d = 8.8E-04;
Te = 633;
b = 2.303E+03;
taw = 1 - (drod/l);
e = 0.73;
Pg = 3;
S = (Drod/(2*d))*Pg;
P = -1.2*log10(S)+23.31;
y = e/(2-e);
u = 1 + (((1-taw))*((1/y)+(1/taw)-1));
c = 0.5*(1 + u - (sqrt((u -1)*(u+3))));
f = (1+c)*((e/(2-e))+(taw*(1-taw)/(1-taw*c)))*(l/drod);
F = 4*f*sigma*drod;
K = (Nrod/pi)*((1/F)+(0.5*(((3-e)/e)*(1/(r0*sigma)))));
Q= 307.59 - 190.96*(log(t/24)).^0.24;
Tc = ((Te^4)+K.*Q).^0.25;
disp(Tc)
tr = t0*exp((b*P)/Tc);
dCDF = 1/tr;
end
回答（1 个）
e = [0.4,0.5,0.6,0.7,0.8,0.9];
Pg = [3,4,5];
tistart = 500;
tEnd = 7200;
CDF_in = 0;
for i = 1:numel(e)
for j = 1:numel(Pg)
[sol(i,j),fval(i,j)] = fsolve(@(x)fun(x,tEnd,CDF_in,e(i),Pg(j)),tistart,optimset('TolX',1e-16,'TolFun',1e-16,'Display','none'));
[tsol{i,j},CDFsol{i,j}] = ode45(@(t,CDF)firstODEfun(t,CDF,e(i),Pg(j)),[sol(i,j) tEnd], CDF_in );
end
end
hold on
plot(tsol{1,1},CDFsol{1,1}(:,1))
plot(tsol{3,2},CDFsol{3,2}(:,1))
plot(tsol{end,end},CDFsol{end,end}(:,1))
hold off
grid on
function res = fun(ti,tEnd,CDF_in,e,Pg)
[tsol, CDFsol] = ode45(@(t,CDF) firstODEfun(t,CDF,e,Pg), [ti tEnd], CDF_in);
res = CDFsol(end,1) - 0.05;
end
function dCDF = firstODEfun(t,~,e,Pg)
t0 = 2.512E-018;
r0 = 8.275E-02;
Nrod = 91;
l = 1.59E-02;
drod = 0.0135;
Drod = 0.01174;
sigma = 5.67E-08;
d = 8.8E-04;
Te = 633;
b = 2.303E+03;
taw = 1 - (drod/l);
%e = 0.73;
%Pg = 3;
S = (Drod/(2*d))*Pg;
P = -1.2*log10(S)+23.31;
y = e/(2-e);
u = 1 + (((1-taw))*((1/y)+(1/taw)-1));
c = 0.5*(1 + u - (sqrt((u -1)*(u+3))));
f = (1+c)*((e/(2-e))+(taw*(1-taw)/(1-taw*c)))*(l/drod);
F = 4*f*sigma*drod;
K = (Nrod/pi)*((1/F)+(0.5*(((3-e)/e)*(1/(r0*sigma)))));
Q= 307.59 - 190.96*(log(t/24)).^0.24;
Tc = ((Te^4)+K.*Q).^0.25;
tr = t0*exp((b*P)/Tc);
dCDF = 1/tr;
end
3 个评论
Alireza Mofidi
2023-5-17
Torsten
2023-5-17
I adjusted the code above.
Alireza Mofidi
2023-5-17
类别
在 帮助中心 和 File Exchange 中查找有关 Programming 的更多信息
2023-5-17
2023-5-17
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
