Iteration over initial condition in ode45

12 ビュー (過去 30 日間)
Alireza Mofidi
Alireza Mofidi 2023 年 5 月 17 日
コメント済み: Alireza Mofidi 2023 年 5 月 17 日
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 件)

Torsten
Torsten 2023 年 5 月 17 日
編集済み: Torsten 2023 年 5 月 17 日
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 件のコメント
Torsten
Torsten 2023 年 5 月 17 日
I adjusted the code above.
Alireza Mofidi
Alireza Mofidi 2023 年 5 月 17 日
I appreciate your help.

サインインしてコメントする。

カテゴリ

Help Center および File ExchangeOrdinary Differential Equations についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by