ODE system with variable parameter, Boundary conditions and domain

Hi all

is it possible to solve an ODE system with these characteristics? The problem can be summerized in this way, calling:

Domain=D (integration domain it is a certain set of row vectors. Fo example 0:0.1:1; 0:0.1:2; 0:0.1:5 etc.)

Parameter=P (inside one of the ode equations, it is a column vector; its length is equal to Domain number of set: in D1=> P=2, in D2=> P=4 etc.)

Bc=B (column vectors that rapresent boundary conditions) B=B1,B2,...Bn

at the first step i want to solve the ode system for D1 and P1 with B1. At the second step D2, P2 and B1 and so on...after having completed these steps i want to solve for D1,P1 with B2 etc..

Below there is my code

xRange is the Domain

Tt0,Pt0 and M01 are the IC (collected in Y0). Ch is the variable parameter. The unknowns are Pt, Tt and M

L=0.2:0.1:0.4;
Nc=0.00001;
xRange=[0:Nc:L];
P0=2.200000000000000e+05;
T0=300;
M01=0.0001:0.00001:0.5;
for i=1:length(M01)

Tt0(i)=T0*(1+((gamma-1)/2)*M01(i)^2);

Pt0(i)=P0*(1+((gamma-1)/2)*M01(i)^2)^(gamma/(gamma-1));

end
Y0=[Tt0;Pt0;M01];
[xSol,YSol]=ode23(@chambsinglebobb,xRange,Y0);


function dYdx=tubosinglebobb(x,Y)
gamma=1.667;

hradialchamber=3.5;%mm
Hradialchamber=hradialchamber/1000;%m 
Dexttantalum=0.001500000000000;
Aeffettiva=(Dexttantalum*(Hradialchamber))-(pi*(Dexttantalum^2)/4);%m^2
Perimetro=(Dexttantalum+Hradialchamber+2*pi*(Dexttantalum/2));
Tw=803;
fLc=0.006816427132907;
Tt=Y(1);
Pt=Y(2);
M=Y(3);
Ch=[0.0104838041220848,0.00821863791965557,0.00729962973207637];

dTtdx=Ch*(Tt-Tw)*(Perimetro/Aeffettiva);
dPtdx=Pt*((-gamma*((M^2)/2)*fLc*(Perimetro/Aeffettiva))-...
    (gamma*((M^2)/2)*dTtdx*(1/Tt)));
dMdx=M*(((1+((gamma-1)/2)*M^2)/(1-M^2))*((gamma*M^2*fLc*...
    Perimetro)/(2*Aeffettiva))+...
    (0.5*((gamma*M^2)+1)/(1-M^2))*(Ch*(Tt-Tw)*Perimetro*...
    (1+((gamma-1)/2)*M^2))/(Aeffettiva*Tt));

dYdx=[dTtdx;dPtdx;dMdx];
end

Thank you very much for the help

Regards

[EldaEbrithil]

Hi all

for what concern varition of the initial conditions i have solved the problem with this simple code:

Nc=0.00001;
xRange=[0:Nc:0.2];
P0=2.200000e+05;
T0=300;
M01=(0.005);

Tt0=T0;
Pt0=P0;
for i=1:20
[xSol,YSol]=ode23(@chambsinglebobb,xRange,[Tt0;Pt0;M01]);
M01=M01+0.01;
YSol_Tt(:,i)=YSol(:,1);
YSol_Pt(:,i)=YSol(:,2);
YSol_M(:,i)=YSol(:,3);
end
figure(1);set(gcf,'Visible', 'on')
plot(xSol,YSol_Tt)
xlabel('S')
ylabel('Tt ')
figure(2);set(gcf,'Visible', 'on')
plot(xSol,YSol_Pt)
xlabel('S')
ylabel('Pt ')
figure(3);set(gcf,'Visible', 'on')
plot(xSol,YSol_M)
xlabel('S')
ylabel('M ')

Now i have to include domain and parameter variations, but actually i don't know how to do it.

Answer 1

[Bjorn Gustavsson]

This could be done something like this:

Dall = {0:0.1:1,0:0.1:2}; % however many domains you need etc
Y0 = {[Tt01;Pt01;M01],[Tt02;Pt02;M02],[Tt03;Pt03;M03]};
for iDom = 1:numel(Dall)
    xRange = Dall{iDom};
    for iInitial = 1:numel(Y0)
        [xSol{iDom,Iinitial},YSol{iDom,Iinitial}]=ode23(@chambsinglebobb,xRange,Y0{iInitial});
    end
end

HTH

[EldaEbrithil]

Oh this is brilliant!! Thank you very much for the help!! but there is one more problem to solve: what if i want to associate to each domain one value of the parameter inside the ODEs? Pratically the ODEs change at each domain due to the parameter value. Domain1=>Parameter1, Domain2=> Parameter2 etc..

[EldaEbrithil]

I have tried in this way bu with no success

Nc=0.00001;
P0=2.200000e+05;
T0=300;
M01=(0.005);Tt01=T0;Pt01=P0;
M02=(0.006);Tt02=T0;Pt02=P0;
M03=(0.007);Tt03=T0;Pt03=P0;
Dall = {0:Nc:L(1),0:Nc:L(2)}; 
Y0 = {[Tt01;Pt01;M01],[Tt02;Pt02;M02],[Tt03;Pt03;M03]};
Ch_cell={ChL(1),ChL(2)};
for iDom = 1:numel(Dall)
    xRange = Dall{iDom};
    Ch=Ch_cell{iDom};
    for iInitial = 1:numel(Y0)
    [xSol{iDom,iInitial},YSol{iDom,iInitial}]=ode23(@chambsinglebobb,xRange,Y0{iInitial},Ch);
    end
end
function dYdx=tubosinglebobb(x,Y,Ch)
gamma=1.667;
hradialchamber=3.5;%mm
Hradialchamber=hradialchamber/1000;%m 
Dexttantalum=0.001500000000000;
Aeffettiva=(Dexttantalum*(Hradialchamber))-(pi*(Dexttantalum^2)/4);%m^2
Perimetro=(Dexttantalum+Hradialchamber+2*pi*(Dexttantalum/2));
Tw=803;
fLc=0.006816427132907;
Tt=Y(1);
Pt=Y(2);
M=Y(3);

dTtdx=Ch*(Tw-Tt)*(Perimetro/Aeffettiva);
dPtdx=Pt*((-gamma*((M^2)/2)*fLc*(Perimetro/Aeffettiva))-...
    (gamma*((M^2)/2)*dTtdx*(1/Tt)));
dMdx=M*(((1+((gamma-1)/2)*M^2)/(1-M^2))*((gamma*(M^2)*fLc*...
    Perimetro)/(2*Aeffettiva))+...
    ((0.5*((gamma*M^2)+1))*(1+((gamma-1)/2)*M^2)/(1-M^2))*...
    (Ch*(Tw-Tt)*Perimetro)/(Aeffettiva*Tt));

dYdx=[dTtdx;dPtdx;dMdx];
end

[Bjorn Gustavsson]

Well, then you should modify your ODE-function a bit, something like this:

function dYdx=tubosinglebobb(x,Y,Par2vary)
% as before but now with the Par2vary input too..

then simply modify the ode23-call to something like this:

[xSol{iDom,iInitial},YSol{iDom,iInitial}]=ode23(@(t,Y) chambsinglebobb(t,Y,CurrentPar),xRange,Y0{iInitial},Ch);

Where you can vary CurrentPar in every step of the loops.

HTH

[EldaEbrithil]

ok now it works, i have changed a line in this:

for iDom = 1:numel(Dall)
    xRange = Dall{iDom};
    Ch=Ch_cell{iDom};
    for iInitial = 1:numel(Y0)
        [xSol{iDom,iInitial},YSol{iDom,iInitial}]=ode23(@(x,Y)...
            chambsinglebobb(x,Y,Ch) ,xRange,Y0{iInitial});
    end
end

Do you agree with that? Is there a more elegant way to write these three vector cell below?

Dall = {0:Nc:L(1),0:Nc:L(2)};

Y0 = {[Tt01;Pt01;M01],[Tt02;Pt02;M02],[Tt03;Pt03;M03]};

Ch_cell={ChL(1),ChL(2)};

[EldaEbrithil]

Sorry we just answer at the same time xD

[Bjorn Gustavsson]

That's good, now my job is done!

[EldaEbrithil]

Thank you very much!!