Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Using pdepe to solve a system of both pdes and odes

Subject: Using pdepe to solve a system of both pdes and odes

From: Layla

Date: 10 Apr, 2011 04:28:04

Message: 1 of 14

Hello,

I'm trying to solve a system of five 5 differential equations: 2 ODEs and 3 PDEs. Pdepe keeps giving me the error "This DAE appears to be of index greater than 1" and I'm not sure if it's a bug in my code or if the problem I'm trying to solve isn't possible in Matlab. My code is below and any advice or insights you could give on how to proceed would be much appreciated.

Thank you,
Layla

-----
function solve_pde
m = 2;
r = 0:0.05:2;
t = 0:10:200;

sol = pdepe(m,@pde,@pde_ic,@pde_bc,r,t);
u1 = sol(:,:,1);
u2 = sol(:,:,2);
u3 = sol(:,:,3);
u4 = sol(:,:,4);
u5 = sol(:,:,5);
-----
function [c,f,s] = pde(r,t,u,DuDr)

a2 = 4e-5;
a5 = 0.1;
g1 = 50;
d2 = 1;
d3 = 70;
d5 = 100;
g5 = 50;

c = [1; 1; 1; 1; 1];
f = [0; a2*(1-u(1))*(1-u(4))*u(2); u(3); 0; a5*u(5)].*DuDr;
s = [u(1)*(1-u(1))-g1*u(3)*u(1); d2*u(2)*(1-u(2)); d3*(u(2)-u(3)); -u(5)*u(4); d5*u(1)*u(2)-g5*u(5)];

-----
function u0 = pde_ic(r)

if 0<= r<0.04
    u0 = [0.01; 1; 0; 1; 0];
else
    u0 = [1; 0; 0; 1; 0];
end

-----
function [pl,ql,pr,qr] = pde_bc(xl,ul,xr,ur,t)

g1 = 50;
d5 = 100;
g5 = 50;

pl = [0; 0; 0; ul(4); 0];
ql = [1; 0; 0; ul(4); 0];
pr = [1-g1; 1; 1; 0; d5*(1-g1)/g5];
qr = [0; 1; 1; ul(4); 0];

Subject: Using pdepe to solve a system of both pdes and odes

From: Grzegorz Knor

Date: 10 Apr, 2011 08:00:21

Message: 2 of 14

Are you sure of boundary conditions?

If you change them for example in this way:
pl = ul;
ql = [0; 0; 0; 0; 0];
pr = ur;
qr = [0; 0; 0; 0; 0];

it works.
So please check it first.

best regards
Grzegorz

Subject: Using pdepe to solve a system of both pdes and odes

From: Layla

Date: 11 Apr, 2011 04:25:04

Message: 3 of 14

"Grzegorz Knor" wrote in message <inro2l$oq8$1@fred.mathworks.com>...
> Are you sure of boundary conditions?
>
> If you change them for example in this way:
> pl = ul;
> ql = [0; 0; 0; 0; 0];
> pr = ur;
> qr = [0; 0; 0; 0; 0];
>
> it works.
> So please check it first.
>
> best regards
> Grzegorz

Thank you Grzegorz! The boundaries did indeed need changing (as did the PDE statements). Thanks!

Subject: Using pdepe to solve a system of both pdes and odes

From: Ashews

Date: 3 Jun, 2011 09:34:06

Message: 4 of 14

"Layla" wrote in message <intvr0$mem$1@fred.mathworks.com>...
> "Grzegorz Knor" wrote in message <inro2l$oq8$1@fred.mathworks.com>...
> > Are you sure of boundary conditions?
> >
> > If you change them for example in this way:
> > pl = ul;
> > ql = [0; 0; 0; 0; 0];
> > pr = ur;
> > qr = [0; 0; 0; 0; 0];
> >
> > it works.
> > So please check it first.
> >
> > best regards
> > Grzegorz
>
> Thank you Grzegorz! The boundaries did indeed need changing (as did the PDE statements). Thanks!

Hi I am trying to solve a PDE and ODEs simultaneously.

Can you please provide me your actual ODEs and PDEs so that I can see how you have modified it to suite the pdepe fuction handler?

Also if possible can you suggest me know to solve my two equations?

Below is are the equations, I need to solve for Tchill and Tevap.



ODEs:

Xdes_dt = ((15*Dso*exp(-Ea/(R*Tcondenser)))*(Xbar_cond - Xdes_cal)/(Rp*10^2));
 
    Xads_dt = ((15*Dso*exp(-Ea/(R*Tevaporator)))*(Xbar_evap - Xads_cal)/(Rp*10^2));

Tevap_dt = ((-hfg*Msg/Mcp_evap)* Xads_dt))+((hf*Msg/Mcp_evap)*
Xdes_dt)+((UA_chill/Mcp_evap)*(Tchill-Tevap));

PDE:

   Tchill_dt = ((-Uf_chill)*(dTchill_dz))+ ((lamda_chill/(den_chill*Cp_chill))*(d2Tchill_dz2) -(UA_chill/(Vol_chill*den_chill*Cp_chill))*(Tchill-Tevap)


I am really stuck in this. Please help me.

Subject: Using pdepe to solve a system of both pdes and odes

From: Ashews

Date: 3 Jun, 2011 09:35:07

Message: 5 of 14

"Layla" wrote in message <intvr0$mem$1@fred.mathworks.com>...
> "Grzegorz Knor" wrote in message <inro2l$oq8$1@fred.mathworks.com>...
> > Are you sure of boundary conditions?
> >
> > If you change them for example in this way:
> > pl = ul;
> > ql = [0; 0; 0; 0; 0];
> > pr = ur;
> > qr = [0; 0; 0; 0; 0];
> >
> > it works.
> > So please check it first.
> >
> > best regards
> > Grzegorz
>
> Thank you Grzegorz! The boundaries did indeed need changing (as did the PDE statements). Thanks!

Hi I am trying to solve a PDE and ODEs simultaneously.

Can you please provide me your actual ODEs and PDEs so that I can see how you have modified it to suite the pdepe fuction handler?

Also if possible can you suggest me know to solve my two equations?

Below is are the equations, I need to solve for Tchill and Tevap.



ODEs:

Xdes_dt = ((15*Dso*exp(-Ea/(R*Tcondenser)))*(Xbar_cond - Xdes_cal)/(Rp*10^2));
 
    Xads_dt = ((15*Dso*exp(-Ea/(R*Tevaporator)))*(Xbar_evap - Xads_cal)/(Rp*10^2));

Tevap_dt = ((-hfg*Msg/Mcp_evap)* Xads_dt))+((hf*Msg/Mcp_evap)*
Xdes_dt)+((UA_chill/Mcp_evap)*(Tchill-Tevap));

PDE:

   Tchill_dt = ((-Uf_chill)*(dTchill_dz))+ ((lamda_chill/(den_chill*Cp_chill))*(d2Tchill_dz2) -(UA_chill/(Vol_chill*den_chill*Cp_chill))*(Tchill-Tevap)


I am really stuck in this. Please help me.

Subject: Using pdepe to solve a system of both pdes and odes

From: Ashews

Date: 3 Jun, 2011 09:37:02

Message: 6 of 14

"Layla" wrote in message <intvr0$mem$1@fred.mathworks.com>...
> "Grzegorz Knor" wrote in message <inro2l$oq8$1@fred.mathworks.com>...
> > Are you sure of boundary conditions?
> >
> > If you change them for example in this way:
> > pl = ul;
> > ql = [0; 0; 0; 0; 0];
> > pr = ur;
> > qr = [0; 0; 0; 0; 0];
> >
> > it works.
> > So please check it first.
> >
> > best regards
> > Grzegorz
>
> Thank you Grzegorz! The boundaries did indeed need changing (as did the PDE statements). Thanks!

Hi I am trying to solve a PDE and ODEs simultaneously.

Can you please provide me your actual ODEs and PDEs so that I can see how you have modified it to suite the pdepe fuction handler?

Also if possible can you suggest me know to solve my two equations?

Below is are the equations, I need to solve for Tchill and Tevap.



ODEs:

Xdes_dt = ((15*Dso*exp(-Ea/(R*Tcondenser)))*(Xbar_cond - Xdes_cal)/(Rp*10^2));
 
    Xads_dt = ((15*Dso*exp(-Ea/(R*Tevaporator)))*(Xbar_evap - Xads_cal)/(Rp*10^2));

Tevap_dt = ((-hfg*Msg/Mcp_evap)* Xads_dt))+((hf*Msg/Mcp_evap)*
Xdes_dt)+((UA_chill/Mcp_evap)*(Tchill-Tevap));

PDE:

   Tchill_dt = ((-Uf_chill)*(dTchill_dz))+ ((lamda_chill/(den_chill*Cp_chill))*(d2Tchill_dz2) -(UA_chill/(Vol_chill*den_chill*Cp_chill))*(Tchill-Tevap)


I am really stuck in this. Please help me.

Subject: Using pdepe to solve a system of both pdes and odes

From: Ashews

Date: 3 Jun, 2011 09:45:05

Message: 7 of 14

"Layla" wrote in message <intvr0$mem$1@fred.mathworks.com>...
> "Grzegorz Knor" wrote in message <inro2l$oq8$1@fred.mathworks.com>...
> > Are you sure of boundary conditions?
> >
> > If you change them for example in this way:
> > pl = ul;
> > ql = [0; 0; 0; 0; 0];
> > pr = ur;
> > qr = [0; 0; 0; 0; 0];
> >
> > it works.
> > So please check it first.
> >
> > best regards
> > Grzegorz
>
> Thank you Grzegorz! The boundaries did indeed need changing (as did the PDE statements). Thanks!



Hi I am trying to solve a PDE and ODEs simultaneously.

Can you please provide me your actual ODEs and PDEs so that I can see how you have modified it to suite the pdepe fuction handler?

Also if possible can you suggest me know to solve my two equations?

Below is are the equations, I need to solve for Tchill and Tevap.



ODEs:

Xdes_dt = ((15*Dso*exp(-Ea/(R*Tcondenser)))*(Xbar_cond - Xdes_cal)/(Rp*10^2));
 
    Xads_dt = ((15*Dso*exp(-Ea/(R*Tevaporator)))*(Xbar_evap - Xads_cal)/(Rp*10^2));

Tevap_dt = ((-hfg*Msg/Mcp_evap)* Xads_dt))+((hf*Msg/Mcp_evap)*
Xdes_dt)+((UA_chill/Mcp_evap)*(Tchill-Tevap));

PDE:

   Tchill_dt = ((-Uf_chill)*(dTchill_dz))+ ((lamda_chill/(den_chill*Cp_chill))*(d2Tchill_dz2) -(UA_chill/(Vol_chill*den_chill*Cp_chill))*(Tchill-Tevap)


I am really stuck in this. Please help me.

Subject: Using pdepe to solve a system of both pdes and odes

From: Ashews

Date: 3 Jun, 2011 09:58:03

Message: 8 of 14

"Ashews " <mechashews@yahoo.co.in> wrote in message <isaaf1$e2b$1@newscl01ah.mathworks.com>...
> "Layla" wrote in message <intvr0$mem$1@fred.mathworks.com>...
> > "Grzegorz Knor" wrote in message <inro2l$oq8$1@fred.mathworks.com>...
> > > Are you sure of boundary conditions?
> > >
> > > If you change them for example in this way:
> > > pl = ul;
> > > ql = [0; 0; 0; 0; 0];
> > > pr = ur;
> > > qr = [0; 0; 0; 0; 0];
> > >
> > > it works.
> > > So please check it first.
> > >
> > > best regards
> > > Grzegorz
> >
> > Thank you Grzegorz! The boundaries did indeed need changing (as did the PDE statements). Thanks!
>
>
>
> Hi I am trying to solve a PDE and ODEs simultaneously.
>
> Can you please provide me your actual ODEs and PDEs so that I can see how you have modified it to suite the pdepe fuction handler?
>
> Also if possible can you suggest me know to solve my two equations?
>
> Below is are the equations, I need to solve for Tchill and Tevap.
>
>
>
> ODEs:
>
> Xdes_dt = ((15*Dso*exp(-Ea/(R*Tcondenser)))*(Xbar_cond - Xdes_cal)/(Rp*10^2));
>
> Xads_dt = ((15*Dso*exp(-Ea/(R*Tevaporator)))*(Xbar_evap - Xads_cal)/(Rp*10^2));
>
> Tevap_dt = ((-hfg*Msg/Mcp_evap)* Xads_dt))+((hf*Msg/Mcp_evap)*
> Xdes_dt)+((UA_chill/Mcp_evap)*(Tchill-Tevap));
>
> PDE:
>
> Tchill_dt = ((-Uf_chill)*(dTchill_dz))+ ((lamda_chill/(den_chill*Cp_chill))*(d2Tchill_dz2) -(UA_chill/(Vol_chill*den_chill*Cp_chill))*(Tchill-Tevap)
>
>
> I am really stuck in this. Please help me.

I missed the boundary conditions :


Intial Conditions @ time t = 0:
Xdes = Xads = 0
Tchill = Tevap = 303

Boundary Conditions for PDE:

Tchill(z = 0,t) = 287.8

dTchill/dz(z = z1,t) = 0

Subject: Using pdepe to solve a system of both pdes and odes

From: salman

Date: 3 Jun, 2011 10:18:05

Message: 9 of 14

"Ashews " <mechashews@yahoo.co.in> wrote in message <isab7b$fre$1@newscl01ah.mathworks.com>...
> "Ashews " <mechashews@yahoo.co.in> wrote in message <isaaf1$e2b$1@newscl01ah.mathworks.com>...
> > "Layla" wrote in message <intvr0$mem$1@fred.mathworks.com>...
> > > "Grzegorz Knor" wrote in message <inro2l$oq8$1@fred.mathworks.com>...
> > > > Are you sure of boundary conditions?
> > > >
> > > > If you change them for example in this way:
> > > > pl = ul;
> > > > ql = [0; 0; 0; 0; 0];
> > > > pr = ur;
> > > > qr = [0; 0; 0; 0; 0];
> > > >
> > > > it works.
> > > > So please check it first.
> > > >
> > > > best regards
> > > > Grzegorz
> > >
> > > Thank you Grzegorz! The boundaries did indeed need changing (as did the PDE statements). Thanks!
> >
> >
> >
> > Hi I am trying to solve a PDE and ODEs simultaneously.
> >
> > Can you please provide me your actual ODEs and PDEs so that I can see how you have modified it to suite the pdepe fuction handler?
> >
> > Also if possible can you suggest me know to solve my two equations?
> >
> > Below is are the equations, I need to solve for Tchill and Tevap.
> >
> >
> >
> > ODEs:
> >
> > Xdes_dt = ((15*Dso*exp(-Ea/(R*Tcondenser)))*(Xbar_cond - Xdes_cal)/(Rp*10^2));
> >
> > Xads_dt = ((15*Dso*exp(-Ea/(R*Tevaporator)))*(Xbar_evap - Xads_cal)/(Rp*10^2));
> >
> > Tevap_dt = ((-hfg*Msg/Mcp_evap)* Xads_dt))+((hf*Msg/Mcp_evap)*
> > Xdes_dt)+((UA_chill/Mcp_evap)*(Tchill-Tevap));
> >
> > PDE:
> >
> > Tchill_dt = ((-Uf_chill)*(dTchill_dz))+ ((lamda_chill/(den_chill*Cp_chill))*(d2Tchill_dz2) -(UA_chill/(Vol_chill*den_chill*Cp_chill))*(Tchill-Tevap)
> >
> >
> > I am really stuck in this. Please help me.
>
> I missed the boundary conditions :
>
>
> Intial Conditions @ time t = 0:
> Xdes = Xads = 0
> Tchill = Tevap = 303
>
> Boundary Conditions for PDE:
>
> Tchill(z = 0,t) = 287.8
>
> dTchill/dz(z = z1,t) = 0


Ashews use the MOL approach for solving this. a pdepe can never be used to solve the system containing both PDEs and ODEs. secondly it depends on your PDE, if it is parabolic or elliptic, you can solve it very easily using pdepe, but if it is hyperbolic, like a one dimensional scalar wave equation, then dont use pdepe, use instead the Method of Lines approach (MOL) which discretizes the space derivative along the grid while keeps he time derivative continuous.

Subject: Using pdepe to solve a system of both pdes and odes

From: Ashews

Date: 3 Jun, 2011 10:38:04

Message: 10 of 14

"salman " <salmanabdullah9@gmail.com> wrote in message <isacct$ij3$1@newscl01ah.mathworks.com>...
> "Ashews " <mechashews@yahoo.co.in> wrote in message <isab7b$fre$1@newscl01ah.mathworks.com>...
> > "Ashews " <mechashews@yahoo.co.in> wrote in message <isaaf1$e2b$1@newscl01ah.mathworks.com>...
> > > "Layla" wrote in message <intvr0$mem$1@fred.mathworks.com>...
> > > > "Grzegorz Knor" wrote in message <inro2l$oq8$1@fred.mathworks.com>...
> > > > > Are you sure of boundary conditions?
> > > > >
> > > > > If you change them for example in this way:
> > > > > pl = ul;
> > > > > ql = [0; 0; 0; 0; 0];
> > > > > pr = ur;
> > > > > qr = [0; 0; 0; 0; 0];
> > > > >
> > > > > it works.
> > > > > So please check it first.
> > > > >
> > > > > best regards
> > > > > Grzegorz
> > > >
> > > > Thank you Grzegorz! The boundaries did indeed need changing (as did the PDE statements). Thanks!
> > >
> > >
> > >
> > > Hi I am trying to solve a PDE and ODEs simultaneously.
> > >
> > > Can you please provide me your actual ODEs and PDEs so that I can see how you have modified it to suite the pdepe fuction handler?
> > >
> > > Also if possible can you suggest me know to solve my two equations?
> > >
> > > Below is are the equations, I need to solve for Tchill and Tevap.
> > >
> > >
> > >
> > > ODEs:
> > >
> > > Xdes_dt = ((15*Dso*exp(-Ea/(R*Tcondenser)))*(Xbar_cond - Xdes_cal)/(Rp*10^2));
> > >
> > > Xads_dt = ((15*Dso*exp(-Ea/(R*Tevaporator)))*(Xbar_evap - Xads_cal)/(Rp*10^2));
> > >
> > > Tevap_dt = ((-hfg*Msg/Mcp_evap)* Xads_dt))+((hf*Msg/Mcp_evap)*
> > > Xdes_dt)+((UA_chill/Mcp_evap)*(Tchill-Tevap));
> > >
> > > PDE:
> > >
> > > Tchill_dt = ((-Uf_chill)*(dTchill_dz))+ ((lamda_chill/(den_chill*Cp_chill))*(d2Tchill_dz2) -(UA_chill/(Vol_chill*den_chill*Cp_chill))*(Tchill-Tevap)
> > >
> > >
> > > I am really stuck in this. Please help me.
> >
> > I missed the boundary conditions :
> >
> >
> > Intial Conditions @ time t = 0:
> > Xdes = Xads = 0
> > Tchill = Tevap = 303
> >
> > Boundary Conditions for PDE:
> >
> > Tchill(z = 0,t) = 287.8
> >
> > dTchill/dz(z = z1,t) = 0
>
>
> Ashews use the MOL approach for solving this. a pdepe can never be used to solve the system containing both PDEs and ODEs. secondly it depends on your PDE, if it is parabolic or elliptic, you can solve it very easily using pdepe, but if it is hyperbolic, like a one dimensional scalar wave equation, then dont use pdepe, use instead the Method of Lines approach (MOL) which discretizes the space derivative along the grid while keeps he time derivative continuous.

Oh is it so. I read somewhere that we can set the flux term to zero while defining ODE. I had the doubt about the boundary conditions that i need to give for ODE wheen trying to solve using the pdepe.

I had also tried the method of lines approach but I was not successful. Can you check if I need to make any changes to the code, below is the discretisation part

for i = 1:n
    den_chill(i) = refpropm('D','T',Tcondenser,'Q',1,'water');
    Cp_chill(i) = refpropm('C','T',Tcondenser,'Q',1,'water');
    lamda_chill(i) = refpropm('L','T',Tcondenser,'Q',1,'water');
    hf(i) = refpropm('H','T',Tevaporator,'Q',1,'water');
    Xdes_cal_dt(i) = ((15*Dso*exp(-Ea/(R*Tcondenser)))*(Xbar_cond - ...
        Xdes_cal(i))/(Rp*10^2));
 
    Xads_cal_dt(i) = ((15*Dso*exp(-Ea/(R*Tevaporator)))*(Xbar_evap - ...
        Xads_cal(i))/(Rp*10^2));
  

Tevap_dt(i) = ((-hfg*Msg/Mcp_evap)* ((15*Dso*exp(-Ea/(R*Tevaporator)))...
    *(Xbar_evap - Xads_cal(i))/(Rp*10^2)))+((hf(i)*Msg/Mcp_evap)*...
     ((15*Dso*exp(-Ea/(R*Tcondenser)))*(Xbar_cond - ...
       Xdes_cal(i))/(Rp*10^2)))+((UA_chill/Mcp_evap)*(Tchill(i)-Tevap(i)));
    

     if i == 1
         Tchill_dt(i)= 0;
         
     elseif i>1 && i<n
         
     
     Tchill_dt(i) = ((-Uf_chill)*((Tchill(i+1)-Tchill(i-1))/(2*deltaz)))+...
     ((lamda_chill(i)/(den_chill(i)*Cp_chill(i)))*((Tchill(i+1)+Tchill(i-1)...
        -2*Tchill(i))/(deltaz^2)))-(UA_chill/(Vol_chill*den_chill(i)...
          *Cp_chill(i)))*(Tchill(i)-Tevap(i));
      
     elseif i == n
         
     Tchill_dt(i) = ((-2)*(lamda_chill(i)/(den_chill(i)*Cp_chill(i)))...
       *((Tchill(i)+Tchill(i-1))/(deltaz^2)))-...
      (UA_chill/(Vol_chill*den_chill(i)*Cp_chill(i)))*(Tchill(i)-Tevap(i));
             %dwdt(i) = Im*(w(i)-w(i-1))/deltaz^2 - U*(w(i)-v(i));
     end
end

If you want I can send the whole program.

Subject: Using pdepe to solve a system of both pdes and odes

From: salman

Date: 3 Jun, 2011 14:12:19

Message: 11 of 14

"Ashews " <mechashews@yahoo.co.in> wrote in message <isadic$lcc$1@newscl01ah.mathworks.com>...
> "salman " <salmanabdullah9@gmail.com> wrote in message <isacct$ij3$1@newscl01ah.mathworks.com>...
> > "Ashews " <mechashews@yahoo.co.in> wrote in message <isab7b$fre$1@newscl01ah.mathworks.com>...
> > > "Ashews " <mechashews@yahoo.co.in> wrote in message <isaaf1$e2b$1@newscl01ah.mathworks.com>...
> > > > "Layla" wrote in message <intvr0$mem$1@fred.mathworks.com>...
> > > > > "Grzegorz Knor" wrote in message <inro2l$oq8$1@fred.mathworks.com>...
> > > > > > Are you sure of boundary conditions?
> > > > > >
> > > > > > If you change them for example in this way:
> > > > > > pl = ul;
> > > > > > ql = [0; 0; 0; 0; 0];
> > > > > > pr = ur;
> > > > > > qr = [0; 0; 0; 0; 0];
> > > > > >
> > > > > > it works.
> > > > > > So please check it first.
> > > > > >
> > > > > > best regards
> > > > > > Grzegorz
> > > > >
> > > > > Thank you Grzegorz! The boundaries did indeed need changing (as did the PDE statements). Thanks!
> > > >
> > > >
> > > >
> > > > Hi I am trying to solve a PDE and ODEs simultaneously.
> > > >
> > > > Can you please provide me your actual ODEs and PDEs so that I can see how you have modified it to suite the pdepe fuction handler?
> > > >
> > > > Also if possible can you suggest me know to solve my two equations?
> > > >
> > > > Below is are the equations, I need to solve for Tchill and Tevap.
> > > >
> > > >
> > > >
> > > > ODEs:
> > > >
> > > > Xdes_dt = ((15*Dso*exp(-Ea/(R*Tcondenser)))*(Xbar_cond - Xdes_cal)/(Rp*10^2));
> > > >
> > > > Xads_dt = ((15*Dso*exp(-Ea/(R*Tevaporator)))*(Xbar_evap - Xads_cal)/(Rp*10^2));
> > > >
> > > > Tevap_dt = ((-hfg*Msg/Mcp_evap)* Xads_dt))+((hf*Msg/Mcp_evap)*
> > > > Xdes_dt)+((UA_chill/Mcp_evap)*(Tchill-Tevap));
> > > >
> > > > PDE:
> > > >
> > > > Tchill_dt = ((-Uf_chill)*(dTchill_dz))+ ((lamda_chill/(den_chill*Cp_chill))*(d2Tchill_dz2) -(UA_chill/(Vol_chill*den_chill*Cp_chill))*(Tchill-Tevap)
> > > >
> > > >
> > > > I am really stuck in this. Please help me.
> > >
> > > I missed the boundary conditions :
> > >
> > >
> > > Intial Conditions @ time t = 0:
> > > Xdes = Xads = 0
> > > Tchill = Tevap = 303
> > >
> > > Boundary Conditions for PDE:
> > >
> > > Tchill(z = 0,t) = 287.8
> > >
> > > dTchill/dz(z = z1,t) = 0
> >
> >
> > Ashews use the MOL approach for solving this. a pdepe can never be used to solve the system containing both PDEs and ODEs. secondly it depends on your PDE, if it is parabolic or elliptic, you can solve it very easily using pdepe, but if it is hyperbolic, like a one dimensional scalar wave equation, then dont use pdepe, use instead the Method of Lines approach (MOL) which discretizes the space derivative along the grid while keeps he time derivative continuous.
>
> Oh is it so. I read somewhere that we can set the flux term to zero while defining ODE. I had the doubt about the boundary conditions that i need to give for ODE wheen trying to solve using the pdepe.
>
> I had also tried the method of lines approach but I was not successful. Can you check if I need to make any changes to the code, below is the discretisation part
>
> for i = 1:n
> den_chill(i) = refpropm('D','T',Tcondenser,'Q',1,'water');
> Cp_chill(i) = refpropm('C','T',Tcondenser,'Q',1,'water');
> lamda_chill(i) = refpropm('L','T',Tcondenser,'Q',1,'water');
> hf(i) = refpropm('H','T',Tevaporator,'Q',1,'water');
> Xdes_cal_dt(i) = ((15*Dso*exp(-Ea/(R*Tcondenser)))*(Xbar_cond - ...
> Xdes_cal(i))/(Rp*10^2));
>
> Xads_cal_dt(i) = ((15*Dso*exp(-Ea/(R*Tevaporator)))*(Xbar_evap - ...
> Xads_cal(i))/(Rp*10^2));
>
>
> Tevap_dt(i) = ((-hfg*Msg/Mcp_evap)* ((15*Dso*exp(-Ea/(R*Tevaporator)))...
> *(Xbar_evap - Xads_cal(i))/(Rp*10^2)))+((hf(i)*Msg/Mcp_evap)*...
> ((15*Dso*exp(-Ea/(R*Tcondenser)))*(Xbar_cond - ...
> Xdes_cal(i))/(Rp*10^2)))+((UA_chill/Mcp_evap)*(Tchill(i)-Tevap(i)));
>
>
> if i == 1
> Tchill_dt(i)= 0;
>
> elseif i>1 && i<n
>
>
> Tchill_dt(i) = ((-Uf_chill)*((Tchill(i+1)-Tchill(i-1))/(2*deltaz)))+...
> ((lamda_chill(i)/(den_chill(i)*Cp_chill(i)))*((Tchill(i+1)+Tchill(i-1)...
> -2*Tchill(i))/(deltaz^2)))-(UA_chill/(Vol_chill*den_chill(i)...
> *Cp_chill(i)))*(Tchill(i)-Tevap(i));
>
> elseif i == n
>
> Tchill_dt(i) = ((-2)*(lamda_chill(i)/(den_chill(i)*Cp_chill(i)))...
> *((Tchill(i)+Tchill(i-1))/(deltaz^2)))-...
> (UA_chill/(Vol_chill*den_chill(i)*Cp_chill(i)))*(Tchill(i)-Tevap(i));
> %dwdt(i) = Im*(w(i)-w(i-1))/deltaz^2 - U*(w(i)-v(i));
> end
> end
>
> If you want I can send the whole program.

Yes thats right. you need to define the correct boundary or intiail conditions for the MOL approach. if yo uneed to define the intial condition then include the initial cohditiosn in the s vector that is passed to the odesolver(@fname,t,s). further more, you need to discretize the system and solve the ODEs +PDE on every grid point. so for every grid point you have to pass the value of the initial conditions, e.g for a gaussian beam E(z,0)=exp(-z.^2/b.^2); where b is widht of the gaussian pulse. so define this for every grid point. yo uunderstand now?

Subject: Using pdepe to solve a system of both pdes and odes

From: salman

Date: 3 Jun, 2011 14:16:05

Message: 12 of 14

"Ashews " <mechashews@yahoo.co.in> wrote in message <isadic$lcc$1@newscl01ah.mathworks.com>...
> "salman " <salmanabdullah9@gmail.com> wrote in message <isacct$ij3$1@newscl01ah.mathworks.com>...
> > "Ashews " <mechashews@yahoo.co.in> wrote in message <isab7b$fre$1@newscl01ah.mathworks.com>...
> > > "Ashews " <mechashews@yahoo.co.in> wrote in message <isaaf1$e2b$1@newscl01ah.mathworks.com>...
> > > > "Layla" wrote in message <intvr0$mem$1@fred.mathworks.com>...
> > > > > "Grzegorz Knor" wrote in message <inro2l$oq8$1@fred.mathworks.com>...
> > > > > > Are you sure of boundary conditions?
> > > > > >
> > > > > > If you change them for example in this way:
> > > > > > pl = ul;
> > > > > > ql = [0; 0; 0; 0; 0];
> > > > > > pr = ur;
> > > > > > qr = [0; 0; 0; 0; 0];
> > > > > >
> > > > > > it works.
> > > > > > So please check it first.
> > > > > >
> > > > > > best regards
> > > > > > Grzegorz
> > > > >
> > > > > Thank you Grzegorz! The boundaries did indeed need changing (as did the PDE statements). Thanks!
> > > >
> > > >
> > > >
> > > > Hi I am trying to solve a PDE and ODEs simultaneously.
> > > >
> > > > Can you please provide me your actual ODEs and PDEs so that I can see how you have modified it to suite the pdepe fuction handler?
> > > >
> > > > Also if possible can you suggest me know to solve my two equations?
> > > >
> > > > Below is are the equations, I need to solve for Tchill and Tevap.
> > > >
> > > >
> > > >
> > > > ODEs:
> > > >
> > > > Xdes_dt = ((15*Dso*exp(-Ea/(R*Tcondenser)))*(Xbar_cond - Xdes_cal)/(Rp*10^2));
> > > >
> > > > Xads_dt = ((15*Dso*exp(-Ea/(R*Tevaporator)))*(Xbar_evap - Xads_cal)/(Rp*10^2));
> > > >
> > > > Tevap_dt = ((-hfg*Msg/Mcp_evap)* Xads_dt))+((hf*Msg/Mcp_evap)*
> > > > Xdes_dt)+((UA_chill/Mcp_evap)*(Tchill-Tevap));
> > > >
> > > > PDE:
> > > >
> > > > Tchill_dt = ((-Uf_chill)*(dTchill_dz))+ ((lamda_chill/(den_chill*Cp_chill))*(d2Tchill_dz2) -(UA_chill/(Vol_chill*den_chill*Cp_chill))*(Tchill-Tevap)
> > > >
> > > >
> > > > I am really stuck in this. Please help me.
> > >
> > > I missed the boundary conditions :
> > >
> > >
> > > Intial Conditions @ time t = 0:
> > > Xdes = Xads = 0
> > > Tchill = Tevap = 303
> > >
> > > Boundary Conditions for PDE:
> > >
> > > Tchill(z = 0,t) = 287.8
> > >
> > > dTchill/dz(z = z1,t) = 0
> >
> >
> > Ashews use the MOL approach for solving this. a pdepe can never be used to solve the system containing both PDEs and ODEs. secondly it depends on your PDE, if it is parabolic or elliptic, you can solve it very easily using pdepe, but if it is hyperbolic, like a one dimensional scalar wave equation, then dont use pdepe, use instead the Method of Lines approach (MOL) which discretizes the space derivative along the grid while keeps he time derivative continuous.
>
> Oh is it so. I read somewhere that we can set the flux term to zero while defining ODE. I had the doubt about the boundary conditions that i need to give for ODE wheen trying to solve using the pdepe.
>
> I had also tried the method of lines approach but I was not successful. Can you check if I need to make any changes to the code, below is the discretisation part
>
> for i = 1:n
> den_chill(i) = refpropm('D','T',Tcondenser,'Q',1,'water');
> Cp_chill(i) = refpropm('C','T',Tcondenser,'Q',1,'water');
> lamda_chill(i) = refpropm('L','T',Tcondenser,'Q',1,'water');
> hf(i) = refpropm('H','T',Tevaporator,'Q',1,'water');
> Xdes_cal_dt(i) = ((15*Dso*exp(-Ea/(R*Tcondenser)))*(Xbar_cond - ...
> Xdes_cal(i))/(Rp*10^2));
>
> Xads_cal_dt(i) = ((15*Dso*exp(-Ea/(R*Tevaporator)))*(Xbar_evap - ...
> Xads_cal(i))/(Rp*10^2));
>
>
> Tevap_dt(i) = ((-hfg*Msg/Mcp_evap)* ((15*Dso*exp(-Ea/(R*Tevaporator)))...
> *(Xbar_evap - Xads_cal(i))/(Rp*10^2)))+((hf(i)*Msg/Mcp_evap)*...
> ((15*Dso*exp(-Ea/(R*Tcondenser)))*(Xbar_cond - ...
> Xdes_cal(i))/(Rp*10^2)))+((UA_chill/Mcp_evap)*(Tchill(i)-Tevap(i)));
>
>
> if i == 1
> Tchill_dt(i)= 0;
>
> elseif i>1 && i<n
>
>
> Tchill_dt(i) = ((-Uf_chill)*((Tchill(i+1)-Tchill(i-1))/(2*deltaz)))+...
> ((lamda_chill(i)/(den_chill(i)*Cp_chill(i)))*((Tchill(i+1)+Tchill(i-1)...
> -2*Tchill(i))/(deltaz^2)))-(UA_chill/(Vol_chill*den_chill(i)...
> *Cp_chill(i)))*(Tchill(i)-Tevap(i));
>
> elseif i == n
>
> Tchill_dt(i) = ((-2)*(lamda_chill(i)/(den_chill(i)*Cp_chill(i)))...
> *((Tchill(i)+Tchill(i-1))/(deltaz^2)))-...
> (UA_chill/(Vol_chill*den_chill(i)*Cp_chill(i)))*(Tchill(i)-Tevap(i));
> %dwdt(i) = Im*(w(i)-w(i-1))/deltaz^2 - U*(w(i)-v(i));
> end
> end
>
> If you want I can send the whole program.

Moreover, you mentioned the flux term, remember that elliptic and parabolic pdes yo are talking about. for hyperbolic pdes dont use pdepe function. just apply MOL approach correctly and you will get the correct result. if you are dealing with parabolic or elliptic pdes then pdepe is best tool for that.

Subject: Using pdepe to solve a system of both pdes and odes

From: Nikhil Chandrikapure

Date: 4 Apr, 2013 19:39:09

Message: 13 of 14

"Layla" wrote in message <inrbkk$deo$1@fred.mathworks.com>...
> Hello,
>
> I'm trying to solve a system of five 5 differential equations: 2 ODEs and 3 PDEs. Pdepe keeps giving me the error "This DAE appears to be of index greater than 1" and I'm not sure if it's a bug in my code or if the problem I'm trying to solve isn't possible in Matlab. My code is below and any advice or insights you could give on how to proceed would be much appreciated.
>
> Thank you,
> Layla
>
> -----
> function solve_pde
> m = 2;
> r = 0:0.05:2;
> t = 0:10:200;
>
> sol = pdepe(m,@pde,@pde_ic,@pde_bc,r,t);
> u1 = sol(:,:,1);
> u2 = sol(:,:,2);
> u3 = sol(:,:,3);
> u4 = sol(:,:,4);
> u5 = sol(:,:,5);
> -----
> function [c,f,s] = pde(r,t,u,DuDr)
>
> a2 = 4e-5;
> a5 = 0.1;
> g1 = 50;
> d2 = 1;
> d3 = 70;
> d5 = 100;
> g5 = 50;
>
> c = [1; 1; 1; 1; 1];
> f = [0; a2*(1-u(1))*(1-u(4))*u(2); u(3); 0; a5*u(5)].*DuDr;
> s = [u(1)*(1-u(1))-g1*u(3)*u(1); d2*u(2)*(1-u(2)); d3*(u(2)-u(3)); -u(5)*u(4); d5*u(1)*u(2)-g5*u(5)];
>
> -----
> function u0 = pde_ic(r)
>
> if 0<= r<0.04
> u0 = [0.01; 1; 0; 1; 0];
> else
> u0 = [1; 0; 0; 1; 0];
> end
>
> -----
> function [pl,ql,pr,qr] = pde_bc(xl,ul,xr,ur,t)
>
> g1 = 50;
> d5 = 100;
> g5 = 50;
>
> pl = [0; 0; 0; ul(4); 0];
> ql = [1; 0; 0; ul(4); 0];
> pr = [1-g1; 1; 1; 0; d5*(1-g1)/g5];
> qr = [0; 1; 1; ul(4); 0];

Hey,i have 3 ODEs and 2 PDEs problem,
my code is here, please suggest me why is this not working..

function [u1,u2,u3,u4,u5]=solve_pde()
m = 1;
x = 0:0.1:1;
t = 0:1:20;
k1=2.81; k2=3.21*(10^-3); k3=4.11*(10^-6);

%for t=0:1:20
%phi(t+1)=exp(-((0.0001*t)^(0.4516)));
%end


sol = pdepe(m,@pdefun,@pde_ic,@pde_bc,x,t);
u1 = sol(:,:,1);
u2 = sol(:,:,2);
u3 = sol(:,:,3);
u4 = sol(:,:,4);
u5 = sol(:,:,5);


function [c,f,s]=pdefun(x,t,u,DuDx)
    

%a1=-2*k1/phi(t+1);
%a2=2*k1/(phi(t+1)^2);
%a3=-k2/(phi(t+1));
%a4=-2*k3/(phi(t+1));
a1=-2*k1;
a2=2*k1;
a3=-k2;
a4=-2*k3;
a5=3.42;
a6=2.124;

c=[1;1;1;1/a5;1/a6];
f=[0;0;0;DuDx;DuDx];
s=[a1*(u(1)^2)+k2*u(3)*u(4);
   a1*(u(1)^2)+k2*u(3)*u(4);
   a2*(u(1)^2)+a3*u(3)*u(4)+a4*u(1);
   -(a1/a5)*(u(1)^2)-(k2/a5)*u(3)*u(4);
   (k3/a6)*u(1)];
end

function u0 = pde_ic()
u0=[0.4;0.4;5.6;0;0];
end

function [pl,ql,pr,qr] = pde_bc(xl,ul,xr,ur,t)
pl=[0;0;0;0;0];
ql=[1;1;1;1;1];
pr=[ur(1);ur(2);ur(3);1;1];
qr=[0;0;0;0;0];

end

end

Subject: Using pdepe to solve a system of both pdes and odes

From: Nikhil Chandrikapure

Date: 4 Apr, 2013 19:55:06

Message: 14 of 14

"Layla" wrote in message <inrbkk$deo$1@fred.mathworks.com>...
> Hello,
>
> I'm trying to solve a system of five 5 differential equations: 2 ODEs and 3 PDEs. Pdepe keeps giving me the error "This DAE appears to be of index greater than 1" and I'm not sure if it's a bug in my code or if the problem I'm trying to solve isn't possible in Matlab. My code is below and any advice or insights you could give on how to proceed would be much appreciated.
>
> Thank you,
> Layla
>
> -----
> function solve_pde
> m = 2;
> r = 0:0.05:2;
> t = 0:10:200;
>
> sol = pdepe(m,@pde,@pde_ic,@pde_bc,r,t);
> u1 = sol(:,:,1);
> u2 = sol(:,:,2);
> u3 = sol(:,:,3);
> u4 = sol(:,:,4);
> u5 = sol(:,:,5);
> -----
> function [c,f,s] = pde(r,t,u,DuDr)
>
> a2 = 4e-5;
> a5 = 0.1;
> g1 = 50;
> d2 = 1;
> d3 = 70;
> d5 = 100;
> g5 = 50;
>
> c = [1; 1; 1; 1; 1];
> f = [0; a2*(1-u(1))*(1-u(4))*u(2); u(3); 0; a5*u(5)].*DuDr;
> s = [u(1)*(1-u(1))-g1*u(3)*u(1); d2*u(2)*(1-u(2)); d3*(u(2)-u(3)); -u(5)*u(4); d5*u(1)*u(2)-g5*u(5)];
>
> -----
> function u0 = pde_ic(r)
>
> if 0<= r<0.04
> u0 = [0.01; 1; 0; 1; 0];
> else
> u0 = [1; 0; 0; 1; 0];
> end
>
> -----
> function [pl,ql,pr,qr] = pde_bc(xl,ul,xr,ur,t)
>
> g1 = 50;
> d5 = 100;
> g5 = 50;
>
> pl = [0; 0; 0; ul(4); 0];
> ql = [1; 0; 0; ul(4); 0];
> pr = [1-g1; 1; 1; 0; d5*(1-g1)/g5];
> qr = [0; 1; 1; ul(4); 0];


can you post here your word to word code..please i need it

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us