Solving a system of ODE over different intervals (with different conditions on parameters)
32 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Lenilein
el 14 de En. de 2019
Comentada: Lenilein
el 15 de En. de 2019
Good Evening dear Matlab Community,
Now that my first part of code works, I have the ambition to make it a little complicated but I can't figure out how to implement my plan in matlab:
in the code below you can see that I'm solving a system of two differential equations. Now my challenge is that I would like to solve it with different conditions on certain parameters over the full interval l. I would like to divide the full interval l into sections i each itself divided into two zones, j.
Each section i starts with a conductive zone (j=1) of say l=3 meter for which:
Aap=lcont
Acp=Aap
hap=17
kap=0.018
And ends with a convection zone (j=2) of l=2 meter for which following is true:
Aap=2*lcont
Acp=0
hap=20
kap=0.02
After this, comes the conductive zone of the next section...
I am confused about how to define indices for the sections and zones and whether I need to include loops within the functions corresponding to my ODE in order to solve my system. I have tried (and failed) by introducing following lines within the script of my ODE functions:
for j=1
Aap=lcont
Acp=Aap
hap=17
kap=0.018
end
for j=2
Aap=2*lcont
Acp=0
hap=20
kap=0.02
end
If I define the array J=[1 2], how can I make sure that the same indice is used to solve both ODE for each section?
Ideally I would even like to be able to build-in variation for the parameters of the conductive zone between two different sections (like for section 1 and j=1, hap=17 and for section 2 and j=1, hap=5) but I guess once I understand how to implement it for one case, I'll be able to extend it to further! :)
[l,y]=ode45(@myODE,[0 45],[0.48,30]);
plot(l,y(:,:))
function dy = myODE(l,y)
global u
u=y(1);
global Tp
Tp=y(2);
dTpdl = myODE1(l,u,Tp);
dudl = myODE2(l,u,Tp);
dy = [dudl;dTpdl];
end
function dudl = myODE2(l,u,Tp)
lcont=2.855;
Psatp = 0.13*exp(18-3800/(Tp+227.03));
Phi= 1-exp(-(47*u^2+0.1*Tp*u^1.06));
Ppw= Psatp*Phi;
dudl=-(kap*Aap)*(Ppw/(Tp+273.15)-0.19);
end
function dTpdl=myODE1(l,u,Tp)
lcont=2.855;
hcp0=0.65;
Psatp = 0.13*exp(18-3800/(Tp+227.03));
Phi= 1-exp(-(47*u^2+0.1*Tp*u^1.06));
Ppw= Psatp*Phi;
DHs = 0.1*(u^1.1)*((Tp+275.15)^2)*(1-Phi)/Phi;
DHvap = 1000*(2501-2.3*Tp);
DHevap = DHs+DHvap;
dTpdl=((hcp0+0.955*u)*(140-Tp)*Acp+hap*(100-Tp)*Aap+DHevap*(-kap*Aap)*(Ppw/(Tp+273)-0.19))/(1.4+4.18*u);
end
3 comentarios
Jan
el 15 de En. de 2019
This seems to mean, that you have to split the integration "time" (the physical meaning does not matter) into intervals like in the 1st example in my answer:
tPool = [0,a,a+b,2*a+b, ...];
Now integrate the function in these intervals and use the final value of one integration as initial value to the next one.
Respuesta aceptada
Jan
el 14 de En. de 2019
Editada: Jan
el 14 de En. de 2019
You can integerate over different time intervals by using a loop:
tPool = [0,1,2,4,8];
y0 = 0; % Your initial value
tAll = [];
yAll = [];
for k = 1:numel(tPool) - 1
% Your parameter depending on the time step:
switch k
case 1
P = 17;
case 2
P = 23;
case 3
P = 10992;
otherwise
P = 0;
end
% Call integrator, provide current parameter by anonymous function:
[t,y] = ode45((t,y) @fcn(t,y,P), [tPool(k:k+1)], y0)
tAll = cat(1, tAll, t); % Append new solution to total solution
yAll = cat(1, yAll, y);
y0 = y(end, 1); % Update the initial value
end
If the zones are not predefined by time steps, but by positions, use a while loop until the final time is reached and an event function set by odeset:
t0 = 0;
tEnd = 8;
y0 = 0; % Your initial value
tAll = [];
yAll = [];
while t0 < tEnd
% Your parameter depending on the time step:
if y(1) < 10
P = 17;
eventP = 20
elseif y(1) < 20
P = 23;
eventP = 30;
else
P = 10992;
eventP = 100;
end
% Create new event function, which stops the integration at the
% next step:
opt = odeset('Events', @(t,y) @myEventsFcn(t,y,eventP));
% Call integrator, provide current parameter by anonymous function:
% tEnd is not reached, because the integration is stopped by the
% event function except for the last step:
[t,y] = ode45((t,y) @fcn(t,y,P), [t0, tEnd], y0, opt)
tAll = cat(1, tAll, t); % Append new solution to total solution
yAll = cat(1, yAll, y);
t0 = t(end);
y0 = y(end, 1); % Update the initial value
end
function [position,isterminal,direction] = yourEventFcn(t,y,eventP)
position = y(1) - eventP; % The value that we want to be zero
isterminal = 1; % Halt integration
direction = 0; % The zero can be approached from either direction
end
Both examples are just a demonstration and they do not run!!!
Más respuestas (0)
Ver también
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!