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

Réponses (1)

Torsten
Torsten le 17 Mai 2023
Modifié(e) : Torsten le 17 Mai 2023
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 commentaires

Alireza Mofidi
Alireza Mofidi le 17 Mai 2023
Thank you very much Torsten for your help. Just one more question, How can I get the CDF solutions vs Time. Currently, it doesn't write the CDFsol vs time.
Torsten
Torsten le 17 Mai 2023
I adjusted the code above.
Alireza Mofidi
Alireza Mofidi le 17 Mai 2023
I appreciate your help.

Connectez-vous pour commenter.

Catégories

En savoir plus sur Programming dans Centre d'aide et File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by