Thread Subject: Can I use pdepe( ) function to solve several first-order PDEs?

Hi, all,

   I want to solve several first-order PDEs as following:

    Du1/Dt = m - u1 - u1*u2/(1+r*u2)
         0 = Du2/Dx - u1*u2/((1+r*u2))
         0 = Du3/Dx + a*u1*u2/((1+r*u2))

   Can I use pdepe function to solve it? since pdepe
requires second-order derivative over x (D^2u/Dt^2) must be
included in f().

   If I can use pdepe( ) to do the work. then another
question: I can only get the boundary values of u1,u2 and
u3 at x = 0, where x has a range from 0 to L(e.g. L=100),
but I can not get the values at x = L. How do i deal with
xr, ur...?

   Or any other methods in Matlab can be used to solve
these equations ?

   Any suggestions are greatly appreciated.

   Many thanks.

Jinlong

>
> Hi, all,
>
> I want to solve several first-order PDEs as
> as following:
>
> Du1/Dt = m - u1 - u1*u2/(1+r*u2)
> 0 = Du2/Dx - u1*u2/((1+r*u2))
> 0 = Du3/Dx + a*u1*u2/((1+r*u2))
>
> Can I use pdepe function to solve it? since pdepe
> requires second-order derivative over x (D^2u/Dt^2)
> must be
> included in f().
>

You shouldn't use pdepe to solve this system.
Formulate it as a DAE-system and use ODE15S or
ODE23T instead.
Concerning the formulation:
1. Write down (eq1) for each spatial grid point
x_1,...,.x_n:
du1(x_i)/dt = m - u1(x_i)-u1(x_i)*u2(x_i)/(1+r*u2(xi))
(2<=i<=n)
2. Substitute the spatial derivates in (eq2) and
(eq3) by finite difference approximations:
du2(xi)/dx ~ (u2(x(i+1))-u2(x(i-1))/(x(i+1)-x(i-1))
(2<=i<=n-1)
du2(xn)/dx ~ (u2(xn)-u2(x(n-1))/(xn-x(n-1))
(or the one-sided difference quotient involving
u2(x(n-2)),u2(x(n-1)) and u2(xn) which approximates
du2(xn)/dx with second order)
(same for u3)
and add them as algebraic equations.
This should give you n-1 differential and 2*(n-1)
algebraic equations to be solved.
(I assume u1(x1),u2(x1) and u3(x1) are known).


> If I can use pdepe( ) to do the work. then another
> her
> question: I can only get the boundary values of u1,u2
> and
> u3 at x = 0, where x has a range from 0 to L(e.g.
> L=100),
> but I can not get the values at x = L. How do i deal
> with
> xr, ur...?
>
> Or any other methods in Matlab can be used to
> to solve
> these equations ?
>
> Any suggestions are greatly appreciated.
>
> Many thanks.
>
> Jinlong

Best wishes
Torsten.

Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote in
message
<21612152.1217488323217.JavaMail.jakarta@nitrogen.mathforum.
org>...
> >
> > Hi, all,
> >
> > I want to solve several first-order PDEs as
> > as following:
> >
> > Du1/Dt = m - u1 - u1*u2/(1+r*u2)
> > 0 = Du2/Dx - u1*u2/((1+r*u2))
> > 0 = Du3/Dx + a*u1*u2/((1+r*u2))
> >
> > Can I use pdepe function to solve it? since pdepe
> > requires second-order derivative over x (D^2u/Dt^2)
> > must be
> > included in f().
> >
>
> You shouldn't use pdepe to solve this system.
> Formulate it as a DAE-system and use ODE15S or
> ODE23T instead.
> Concerning the formulation:
> 1. Write down (eq1) for each spatial grid point
> x_1,...,.x_n:
> du1(x_i)/dt = m - u1(x_i)-u1(x_i)*u2(x_i)/(1+r*u2(xi))
> (2<=i<=n)
> 2. Substitute the spatial derivates in (eq2) and
> (eq3) by finite difference approximations:
> du2(xi)/dx ~ (u2(x(i+1))-u2(x(i-1))/(x(i+1)-x(i-1))
> (2<=i<=n-1)
> du2(xn)/dx ~ (u2(xn)-u2(x(n-1))/(xn-x(n-1))
> (or the one-sided difference quotient involving
> u2(x(n-2)),u2(x(n-1)) and u2(xn) which approximates
> du2(xn)/dx with second order)
> (same for u3)
> and add them as algebraic equations.
> This should give you n-1 differential and 2*(n-1)
> algebraic equations to be solved.
> (I assume u1(x1),u2(x1) and u3(x1) are known).
>

  Thanks. I think you are rigth. just two more questions:
1, I think i should start from 1 but not 2, because if i
starts from 2, then the initial value u1(x1) will not be
used at all since there is no u1(x(i-1)) in the equations.
when i=1, du2(xi)/dx ~ (u2(x(i+1))-u2(xi))/(x(i+1)-xi).

2, you said "and add them as algebraic equations." does it
mean just add them into my nested functions? because I
found a similar example on Brusselator system in Matlab help
, in this example, however, all equations in the nested
function are ODE equations. I wonder whether it is allowed
to put ODE and algebraic equations together in the nested
function.
 

>
> > If I can use pdepe( ) to do the work. then another
> > her
> > question: I can only get the boundary values of u1,u2
> > and
> > u3 at x = 0, where x has a range from 0 to L(e.g.
> > L=100),
> > but I can not get the values at x = L. How do i deal
> > with
> > xr, ur...?
> >
> > Or any other methods in Matlab can be used to
> > to solve
> > these equations ?
> >
> > Any suggestions are greatly appreciated.
> >
> > Many thanks.
> >
> > Jinlong
>
> Best wishes
> Torsten.

> Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote
> in
> message
> <21612152.1217488323217.JavaMail.jakarta@nitrogen.math
> forum.
> org>...
> > >
> > > Hi, all,
> > >
> > > I want to solve several first-order PDEs as
> > > as following:
> > >
> > > Du1/Dt = m - u1 - u1*u2/(1+r*u2)
> > > 0 = Du2/Dx - u1*u2/((1+r*u2))
> > > 0 = Du3/Dx + a*u1*u2/((1+r*u2))
> > >
> > > Can I use pdepe function to solve it? since
> pdepe
> > > requires second-order derivative over x
> (D^2u/Dt^2)
> > > must be
> > > included in f().
> > >
> >
> > You shouldn't use pdepe to solve this system.
> > Formulate it as a DAE-system and use ODE15S or
> > ODE23T instead.
> > Concerning the formulation:
> > 1. Write down (eq1) for each spatial grid point
> > x_1,...,.x_n:
> > du1(x_i)/dt = m -
> u1(x_i)-u1(x_i)*u2(x_i)/(1+r*u2(xi))
> > (2<=i<=n)
> > 2. Substitute the spatial derivates in (eq2) and
> > (eq3) by finite difference approximations:
> > du2(xi)/dx ~ (u2(x(i+1))-u2(x(i-1))/(x(i+1)-x(i-1))
> > (2<=i<=n-1)
> > du2(xn)/dx ~ (u2(xn)-u2(x(n-1))/(xn-x(n-1))
> > (or the one-sided difference quotient involving
> > u2(x(n-2)),u2(x(n-1)) and u2(xn) which approximates
> > du2(xn)/dx with second order)
> > (same for u3)
> > and add them as algebraic equations.
> > This should give you n-1 differential and 2*(n-1)
> > algebraic equations to be solved.
> > (I assume u1(x1),u2(x1) and u3(x1) are known).
> >
>
> Thanks. I think you are rigth. just two more
> re questions:
> 1, I think i should start from 1 but not 2, because
> if i
> starts from 2, then the initial value u1(x1) will not
> be
> used at all since there is no u1(x(i-1)) in the
> equations.
> when i=1, du2(xi)/dx ~
> (u2(x(i+1))-u2(xi))/(x(i+1)-xi).

Because you set boundary conditions for u2 and u3 at x=0,
it follows from your equations that the value of
u1(x1) indeed does not influence u2 or u3.
If this is not correct from the physical point of view,
you should reconsider your boundary conditions.
 
>
> 2, you said "and add them as algebraic equations."
> does it
> mean just add them into my nested functions? because
> I
> found a similar example on Brusselator system in
> Matlab help
> , in this example, however, all equations in the
> nested
> function are ODE equations. I wonder whether it is
> allowed
> to put ODE and algebraic equations together in the
> nested
> function.
>

The system equation reads M*y'(t) = f(t,y)

For DAE systems, you have to set and specify the mass
matrix M in which the rows which refer to algebraic
equations are set to zero.

The right-hand side f of the ODE system can be specified
as usual.


>
> >
> > > If I can use pdepe( ) to do the work. then
> another
> > > her
> > > question: I can only get the boundary values of
> u1,u2
> > > and
> > > u3 at x = 0, where x has a range from 0 to L(e.g.
> > > L=100),
> > > but I can not get the values at x = L. How do i
> deal
> > > with
> > > xr, ur...?
> > >
> > > Or any other methods in Matlab can be used to
> > > to solve
> > > these equations ?
> > >
> > > Any suggestions are greatly appreciated.
> > >
> > > Many thanks.
> > >
> > > Jinlong
> >
> > Best wishes
> > Torsten.
>

Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote in
message
<1715863.1217501095914.JavaMail.jakarta@nitrogen.mathforum.o
rg>...
> > Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote
> > in
> > message
> > <21612152.1217488323217.JavaMail.jakarta@nitrogen.math
> > forum.
> > org>...
> > > >
> > > > Hi, all,
> > > >
> > > > I want to solve several first-order PDEs as
> > > > as following:
> > > >
> > > > Du1/Dt = m - u1 - u1*u2/(1+r*u2)
> > > > 0 = Du2/Dx - u1*u2/((1+r*u2))
> > > > 0 = Du3/Dx + a*u1*u2/((1+r*u2))
> > > >
> > > > Can I use pdepe function to solve it? since
> > pdepe
> > > > requires second-order derivative over x
> > (D^2u/Dt^2)
> > > > must be
> > > > included in f().
> > > >
> > >
> > > You shouldn't use pdepe to solve this system.
> > > Formulate it as a DAE-system and use ODE15S or
> > > ODE23T instead.
> > > Concerning the formulation:
> > > 1. Write down (eq1) for each spatial grid point
> > > x_1,...,.x_n:
> > > du1(x_i)/dt = m -
> > u1(x_i)-u1(x_i)*u2(x_i)/(1+r*u2(xi))
> > > (2<=i<=n)
> > > 2. Substitute the spatial derivates in (eq2) and
> > > (eq3) by finite difference approximations:
> > > du2(xi)/dx ~ (u2(x(i+1))-u2(x(i-1))/(x(i+1)-x(i-1))
> > > (2<=i<=n-1)
> > > du2(xn)/dx ~ (u2(xn)-u2(x(n-1))/(xn-x(n-1))
> > > (or the one-sided difference quotient involving
> > > u2(x(n-2)),u2(x(n-1)) and u2(xn) which approximates
> > > du2(xn)/dx with second order)
> > > (same for u3)
> > > and add them as algebraic equations.
> > > This should give you n-1 differential and 2*(n-1)
> > > algebraic equations to be solved.
> > > (I assume u1(x1),u2(x1) and u3(x1) are known).
> > >
> >
> > Thanks. I think you are rigth. just two more
> > re questions:
> > 1, I think i should start from 1 but not 2, because
> > if i
> > starts from 2, then the initial value u1(x1) will not
> > be
> > used at all since there is no u1(x(i-1)) in the
> > equations.
> > when i=1, du2(xi)/dx ~
> > (u2(x(i+1))-u2(xi))/(x(i+1)-xi).
>
> Because you set boundary conditions for u2 and u3 at x=0,
> it follows from your equations that the value of
> u1(x1) indeed does not influence u2 or u3.
> If this is not correct from the physical point of view,
> you should reconsider your boundary conditions.
>
> >
> > 2, you said "and add them as algebraic equations."
> > does it
> > mean just add them into my nested functions? because
> > I
> > found a similar example on Brusselator system in
> > Matlab help
> > , in this example, however, all equations in the
> > nested
> > function are ODE equations. I wonder whether it is
> > allowed
> > to put ODE and algebraic equations together in the
> > nested
> > function.
> >
>
> The system equation reads M*y'(t) = f(t,y)
>
> For DAE systems, you have to set and specify the mass
> matrix M in which the rows which refer to algebraic
> equations are set to zero.
>
> The right-hand side f of the ODE system can be specified
> as usual.

I set up an model for this by using the above method. But
when I run the code, it complains as following:

??? Error using ==> funfun\private\daeic12 at 77
This DAE appears to be of index greater than 1.

what does it mean? can you please give me some guide on how
to solve this problem cause it is really new to me.

Thanks.
>
>
> >
> > >
> > > > If I can use pdepe( ) to do the work. then
> > another
> > > > her
> > > > question: I can only get the boundary values of
> > u1,u2
> > > > and
> > > > u3 at x = 0, where x has a range from 0 to L(e.g.
> > > > L=100),
> > > > but I can not get the values at x = L. How do i
> > deal
> > > > with
> > > > xr, ur...?
> > > >
> > > > Or any other methods in Matlab can be used to
> > > > to solve
> > > > these equations ?
> > > >
> > > > Any suggestions are greatly appreciated.
> > > >
> > > > Many thanks.
> > > >
> > > > Jinlong
> > >
> > > Best wishes
> > > Torsten.
> >

>I set up an model for this by using the above method. >But
>when I run the code, it complains as following:
>
>??? Error using ==> funfun\private\daeic12 at 77
>This DAE appears to be of index greater than 1.
>
>
>what does it mean? can you please give me some guide on >how
>to solve this problem cause it is really new to me.
>
>
>Thanks.


Could you please post your final set of discretized
equations ?

Best wishes
Torsten.

Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote in
message
<18078046.1217592843628.JavaMail.jakarta@nitrogen.mathforum.
org>...
> >I set up an model for this by using the above method.
>But
> >when I run the code, it complains as following:
> >
> >??? Error using ==> funfun\private\daeic12 at 77
> >This DAE appears to be of index greater than 1.
> >
> >
> >what does it mean? can you please give me some guide on
>how
> >to solve this problem cause it is really new to me.
> >
> >
> >Thanks.
>
>
> Could you please post your final set of discretized
> equations ?
>
> Best wishes
> Torsten.

Ok, just paste the key part of the code.

function [g_out,p_out,phi_out,t_out] = soadae1(Opt_in,
I_bias, fs, N)

.....
delta_z = L_amp/N;
M = zeros(3*(N-1),3*(N-1));
   for ii = 1:3:3*(N-1)
       M(ii,ii) = 1;
   end
tspan = 0: 1/fs: (length(Opt_in) - 1)/fs;
y0 = repmat([G_zero; 0; 0],N-1,1);
   
   options = odeset('Mass',M,'absTol',
val_AbsTol,'Vectorized','on','JPattern',jpattern(N));

[t,y] = ode15s(@f,tspan,y0,options);

function out = f(t, y)
tmp = floor(t*fs)+1;
 i = 2;
out(i-1:i+1,:) = [ Coeff1(1) - Coeff1(2)*y(i-1,:) - Coeff1
(3)*y(i-1,:).*y(i,:)./(1 + Gain_comp*y(i,:));...
                      2*delta_z*y(i-1,:).*y(i,:) - (1 +
Gain_comp*y(i,:)).*(y(i+3,:) - P_in(tmp)); ...
                      Linw_en/2*2*delta_z*y(i-1,:) + (1 +
Gain_comp*y(i,:)).*(y(i+4,:) - Phi_in(tmp))];

for i = 5:3:3*(N-1)-4
   out(i-1:i+1,:) = [ Coeff1(1) - Coeff1(2)*y(i-1,:) -
Coeff1(3)*y(i-1,:).*y(i,:)./(1 + Gain_comp*y(i,:));...
                      2*delta_z*y(i-1,:).*y(i,:) - (1 +
Gain_comp*y(i,:)).*(y(i+3,:) - y(i-3,:)); ...
                      Linw_en/2*2*delta_z*y(i-1,:) + (1 +
Gain_comp*y(i,:)).*(y(i+4,:) - y(i-2,:))];

i = 3*(N-1)-1
   out(i-1:i+1,:) = [ Coeff1(1) - Coeff1(2)*y(i-1,:) -
Coeff1(3)*y(i-1,:).*y(i,:)./(1 + Gain_comp*y(i,:));...
                      delta_z*y(i-1,:).*y(i,:) - (1 +
Gain_comp*y(i,:)).*(y(i,:) - y(i-3,:)); ...
                      Linw_en/2*delta_z*y(i-1,:) + (1 +
Gain_comp*y(i,:)).*(y(i+1,:) - y(i-2,:))];
   end %end nested function
end %end function soadae1()

 function S = jpattern(N)
   B = ones(3*(N-1),8);
   B(1:3:3*(N-1)-2,1:2) = zeros(N-1,2);
   B(1:3:3*(N-1)-2,5:8) = zeros(N-1,4);
   B(2:3:3*(N-1)-1,2) = zeros(N-1,1);
   B(2:3:3*(N-1)-1,5:6) = zeros(N-1,2);
   B(2:3:3*(N-1)-1,8) = zeros(N-1,1);
   B(3:3:3*(N-1),1) = zeros(N-1,1);
   B(3:3:3*(N-1),5:7) = zeros(N-1,3);
   S = spdiags(B,-3:4,3*(N-1),3*(N-1));
end


Thank you very much.

>Ok, just paste the key part of the code.
>
>function [g_out,p_out,phi_out,t_out] = soadae1(Opt_in,
>I_bias, fs, N)
>
>
>.....
>delta_z = L_amp/N;
>M = zeros(3*(N-1),3*(N-1));
> for ii = 1:3:3*(N-1)
> M(ii,ii) = 1;
> end
>tspan = 0: 1/fs: (length(Opt_in) - 1)/fs;
>y0 = repmat([G_zero; 0; 0],N-1,1);
>
>
> options = odeset('Mass',M,'absTol',
>val_AbsTol,'Vectorized','on','JPattern',jpattern(N));
>
>
>[t,y] = ode15s(@f,tspan,y0,options);
>
>
>function out = f(t, y)
>tmp = floor(t*fs)+1;
> i = 2;
>out(i-1:i+1,:) = [ Coeff1(1) - Coeff1(2)*y(i-1,:) - >Coeff1
>(3)*y(i-1,:).*y(i,:)./(1 + Gain_comp*y(i,:));...
> 2*delta_z*y(i-1,:).*y(i,:) - (1 +
>Gain_comp*y(i,:)).*(y(i+3,:) - P_in(tmp)); ...
> Linw_en/2*2*delta_z*y(i-1,:) + (1 >+
>Gain_comp*y(i,:)).*(y(i+4,:) - Phi_in(tmp))];
>
>
>for i = 5:3:3*(N-1)-4
> out(i-1:i+1,:) = [ Coeff1(1) - Coeff1(2)*y(i-1,:) -
>Coeff1(3)*y(i-1,:).*y(i,:)./(1 + Gain_comp*y(i,:));...
> 2*delta_z*y(i-1,:).*y(i,:) - (1 +
>Gain_comp*y(i,:)).*(y(i+3,:) - y(i-3,:)); ...
> Linw_en/2*2*delta_z*y(i-1,:) + (1 >+
>Gain_comp*y(i,:)).*(y(i+4,:) - y(i-2,:))];
>
>
>i = 3*(N-1)-1
> out(i-1:i+1,:) = [ Coeff1(1) - Coeff1(2)*y(i-1,:) -
>Coeff1(3)*y(i-1,:).*y(i,:)./(1 + Gain_comp*y(i,:));...
> delta_z*y(i-1,:).*y(i,:) - (1 +
>Gain_comp*y(i,:)).*(y(i,:) - y(i-3,:)); ...
> Linw_en/2*delta_z*y(i-1,:) + (1 +
>Gain_comp*y(i,:)).*(y(i+1,:) - y(i-2,:))];
> end %end nested function
>end %end function soadae1()
>
>
> function S = jpattern(N)
> B = ones(3*(N-1),8);
> B(1:3:3*(N-1)-2,1:2) = zeros(N-1,2);
> B(1:3:3*(N-1)-2,5:8) = zeros(N-1,4);
> B(2:3:3*(N-1)-1,2) = zeros(N-1,1);
> B(2:3:3*(N-1)-1,5:6) = zeros(N-1,2);
> B(2:3:3*(N-1)-1,8) = zeros(N-1,1);
> B(3:3:3*(N-1),1) = zeros(N-1,1);
> B(3:3:3*(N-1),5:7) = zeros(N-1,3);
> S = spdiags(B,-3:4,3*(N-1),3*(N-1));
>end
>
>
>Thank you very much.

I did not check everything in detail, but we overlooked
the fact that our discretization gives no equation
to determine u3(z2).

So maybe for a first test, a first-order approximation
for the spatial derivatives is sufficient:

du1(i)/dt = c1 - c2*u1(i)-c3*u1(i)*u2(i)/(1+g*u2(i))
0 = (u2(i)-u2(i-1))/deltaz - u1(i)*u2(i)/(1+g*u2(i))
0 = (u3(i)-u3(i-1))/deltaz + l/2*u1(i)/(1+g*u2(i))
(2 <= i <= n)

Best wishes
Torsten.

>I did not check everything in detail, but we overlooked
>the fact that our discretization gives no equation
>to determine u3(z2).
>
>So maybe for a first test, a first-order approximation
>for the spatial derivatives is sufficient:
>
>du1(i)/dt = c1 - c2*u1(i)-c3*u1(i)*u2(i)/(1+g*u2(i))
>0 = (u2(i)-u2(i-1))/deltaz - u1(i)*u2(i)/(1+g*u2(i))
>0 = (u3(i)-u3(i-1))/deltaz + l/2*u1(i)/(1+g*u2(i))
>(2 <= i <= n)
>
>Best wishes
>Torsten.

No, I checked it again and now see that there is
(at least formally) a determining equation for each
of the unknowns.

What might be a problem is the equation for u2(n):
(u2(n)-u2(n-1))/deltaz - u1(n)*u2(n)/(1+g*u2(n)) = 0
This results in a quadratic equation in u2(n) which
might not be solvable or not be uniquely solvable.

Because the equations for u1 and u2 are
independent from u3, I would try to solve these two
in a first step separately.

And maybe you can choose the number of grid points
at the beginning such that you can calculate with a full
Jacobian matrix. I often make mistakes when setting
up the pattern of the Jacobian matrix.

Good luck !

Best wishes
Torsten.

Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote in
message
<342330.1217599706924.JavaMail.jakarta@nitrogen.mathforum.or
g>...
> >Ok, just paste the key part of the code.
> >
> >function [g_out,p_out,phi_out,t_out] = soadae1(Opt_in,
> >I_bias, fs, N)
> >
> >
> >.....
> >delta_z = L_amp/N;
> >M = zeros(3*(N-1),3*(N-1));
> > for ii = 1:3:3*(N-1)
> > M(ii,ii) = 1;
> > end
> >tspan = 0: 1/fs: (length(Opt_in) - 1)/fs;
> >y0 = repmat([G_zero; 0; 0],N-1,1);
> >
> >
> > options = odeset('Mass',M,'absTol',
> >val_AbsTol,'Vectorized','on','JPattern',jpattern(N));
> >
> >
> >[t,y] = ode15s(@f,tspan,y0,options);
> >
> >
> >function out = f(t, y)
> >tmp = floor(t*fs)+1;
> > i = 2;
> >out(i-1:i+1,:) = [ Coeff1(1) - Coeff1(2)*y(i-1,:) -
>Coeff1
> >(3)*y(i-1,:).*y(i,:)./(1 + Gain_comp*y(i,:));...
> > 2*delta_z*y(i-1,:).*y(i,:) - (1 +
> >Gain_comp*y(i,:)).*(y(i+3,:) - P_in(tmp)); ...
> > Linw_en/2*2*delta_z*y(i-1,:) + (1
>+
> >Gain_comp*y(i,:)).*(y(i+4,:) - Phi_in(tmp))];
> >
> >
> >for i = 5:3:3*(N-1)-4
> > out(i-1:i+1,:) = [ Coeff1(1) - Coeff1(2)*y(i-1,:) -
> >Coeff1(3)*y(i-1,:).*y(i,:)./(1 + Gain_comp*y(i,:));...
> > 2*delta_z*y(i-1,:).*y(i,:) - (1 +
> >Gain_comp*y(i,:)).*(y(i+3,:) - y(i-3,:)); ...
> > Linw_en/2*2*delta_z*y(i-1,:) + (1
>+
> >Gain_comp*y(i,:)).*(y(i+4,:) - y(i-2,:))];
> >
> >
> >i = 3*(N-1)-1
> > out(i-1:i+1,:) = [ Coeff1(1) - Coeff1(2)*y(i-1,:) -
> >Coeff1(3)*y(i-1,:).*y(i,:)./(1 + Gain_comp*y(i,:));...
> > delta_z*y(i-1,:).*y(i,:) - (1 +
> >Gain_comp*y(i,:)).*(y(i,:) - y(i-3,:)); ...
> > Linw_en/2*delta_z*y(i-1,:) + (1 +
> >Gain_comp*y(i,:)).*(y(i+1,:) - y(i-2,:))];
> > end %end nested function
> >end %end function soadae1()
> >
> >
> > function S = jpattern(N)
> > B = ones(3*(N-1),8);
> > B(1:3:3*(N-1)-2,1:2) = zeros(N-1,2);
> > B(1:3:3*(N-1)-2,5:8) = zeros(N-1,4);
> > B(2:3:3*(N-1)-1,2) = zeros(N-1,1);
> > B(2:3:3*(N-1)-1,5:6) = zeros(N-1,2);
> > B(2:3:3*(N-1)-1,8) = zeros(N-1,1);
> > B(3:3:3*(N-1),1) = zeros(N-1,1);
> > B(3:3:3*(N-1),5:7) = zeros(N-1,3);
> > S = spdiags(B,-3:4,3*(N-1),3*(N-1));
> >end
> >
> >
> >Thank you very much.
>
> I did not check everything in detail, but we overlooked
> the fact that our discretization gives no equation
> to determine u3(z2).
>
> So maybe for a first test, a first-order approximation
> for the spatial derivatives is sufficient:
>
> du1(i)/dt = c1 - c2*u1(i)-c3*u1(i)*u2(i)/(1+g*u2(i))
> 0 = (u2(i)-u2(i-1))/deltaz - u1(i)*u2(i)/(1+g*u2(i))
> 0 = (u3(i)-u3(i-1))/deltaz + l/2*u1(i)/(1+g*u2(i))
> (2 <= i <= n)
>
> Best wishes
> Torsten.

Hi, Torsten,
the new nested function is listed as following:

   function out = f(t, y)
      tmp = floor(t*fs)+1;
   %first-order approximation
   i = 4;
   out(i-3:i-1,:) = [ Coeff1(1) - Coeff1(2)*y(i,:) - Coeff1
(3)*y(i,:).*y(i+1,:)./(1 + Gain_comp*y(i+1,:));...
 delta_z*y(i,:).*y(i+1,:) - (1 + Gain_comp*y(i+1,:)).*(y
(i+1,:) - P_in(tmp)); ...
  Linw_en/2*delta_z*y(i,:) + (1 + Gain_comp*y(i+1,:)).*(y
(i+2,:) - Phi_in(tmp))];

   for i = 7:3:3*(N-1)-2
   out(i-3:i-1,:) = [ Coeff1(1) - Coeff1(2)*y(i,:) - Coeff1
(3)*y(i,:).*y(i+1,:)./(1 + Gain_comp*y(i+1,:));...
   delta_z*y(i,:).*y(i+1,:) - (1 + Gain_comp*y(i+1,:)).*(y
(i+1,:) - y(i-2,:)); ...
   Linw_en/2*delta_z*y(i,:) + (1 + Gain_comp*y(i+1,:)).*(y
(i+2,:) - y(i-1,:))];
   end
end

      The problem now is that output of the nested function
(see varible 'out') is vector of 3(N-2)-by-1,while the size
of y is 3*(N-1)-by-1. Thus it can not be performed by
Matlab.

      Do I need to set the first 3 rows of 'out' to be
zeros? But it seems not correct.

      Please help me have a look.
      Many thanks.

> Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote
> in
> message
> >
> > I did not check everything in detail, but we
> overlooked
> > the fact that our discretization gives no equation
> > to determine u3(z2).
> >
> > So maybe for a first test, a first-order
> approximation
> > for the spatial derivatives is sufficient:
> >
> > du1(i)/dt = c1 -
> c2*u1(i)-c3*u1(i)*u2(i)/(1+g*u2(i))
> > 0 = (u2(i)-u2(i-1))/deltaz -
> u1(i)*u2(i)/(1+g*u2(i))
> > 0 = (u3(i)-u3(i-1))/deltaz + l/2*u1(i)/(1+g*u2(i))
> > (2 <= i <= n)
> >
> > Best wishes
> > Torsten.
>
> Hi, Torsten,
> the new nested function is listed as following:
>
> function out = f(t, y)
> tmp = floor(t*fs)+1;
> %first-order approximation
> i = 4;
> out(i-3:i-1,:) = [ Coeff1(1) - Coeff1(2)*y(i,:) -
> ) - Coeff1
> (3)*y(i,:).*y(i+1,:)./(1 + Gain_comp*y(i+1,:));...
> delta_z*y(i,:).*y(i+1,:) - (1 +
> + Gain_comp*y(i+1,:)).*(y
> (i+1,:) - P_in(tmp)); ...
> Linw_en/2*delta_z*y(i,:) + (1 +
> + Gain_comp*y(i+1,:)).*(y
> (i+2,:) - Phi_in(tmp))];
>
> for i = 7:3:3*(N-1)-2
> out(i-3:i-1,:) = [ Coeff1(1) - Coeff1(2)*y(i,:) -
> ) - Coeff1
> (3)*y(i,:).*y(i+1,:)./(1 + Gain_comp*y(i+1,:));...
> delta_z*y(i,:).*y(i+1,:) - (1 +
> 1 + Gain_comp*y(i+1,:)).*(y
> (i+1,:) - y(i-2,:)); ...
> Linw_en/2*delta_z*y(i,:) + (1 +
> 1 + Gain_comp*y(i+1,:)).*(y
> (i+2,:) - y(i-1,:))];
> end
> end
>
> The problem now is that output of the nested
> nested function
> (see varible 'out') is vector of 3(N-2)-by-1,while
> the size
> of y is 3*(N-1)-by-1. Thus it can not be performed
> by
> Matlab.
>
> Do I need to set the first 3 rows of 'out' to
> ut' to be
> zeros? But it seems not correct.
>
> Please help me have a look.
> Many thanks.
>

It's just a question of enumeration of the variables.
The system

du1(i)/dt = c1 - c2*u1(i)-c3*u1(i)*u2(i)/(1+g*u2(i))
 0 = (u2(i)-u2(i-1))/deltaz -u1(i)*u2(i)/(1+g*u2(i))
 0 = (u3(i)-u3(i-1))/deltaz + l/2*u1(i)/(1+g*u2(i))
(2 <= i <= n)

has 3*(n-1) equations for the 3*(n-1) unknowns
u1(2),...,u1(n),u2(2),...,u2(n),u3(2),...,u3(n).

You have to shift the enumeration for one index
(2 becomes 1, 3 becomes 2 etc.)

Another way is to add the equations
du1(1)/dt = c1 - c2*u1(1)-c3*u1(1)*u2(1)/(1+g*u2(1))
0 = u2(1) - P_in(tmp)
0 = u3(1) - Phi_in(tmp)
and leave the enumeration as it is.

Best wishes
Torsten.

Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote in
message
<15852838.1217677313101.JavaMail.jakarta@nitrogen.mathforum.
org>...
> > Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote
> > in
> > message
> > >
> > > I did not check everything in detail, but we
> > overlooked
> > > the fact that our discretization gives no equation
> > > to determine u3(z2).
> > >
> > > So maybe for a first test, a first-order
> > approximation
> > > for the spatial derivatives is sufficient:
> > >
> > > du1(i)/dt = c1 -
> > c2*u1(i)-c3*u1(i)*u2(i)/(1+g*u2(i))
> > > 0 = (u2(i)-u2(i-1))/deltaz -
> > u1(i)*u2(i)/(1+g*u2(i))
> > > 0 = (u3(i)-u3(i-1))/deltaz + l/2*u1(i)/(1+g*u2(i))
> > > (2 <= i <= n)
> > >
> > > Best wishes
> > > Torsten.
> >
> > Hi, Torsten,
> > the new nested function is listed as following:
> >
> > function out = f(t, y)
> > tmp = floor(t*fs)+1;
> > %first-order approximation
> > i = 4;
> > out(i-3:i-1,:) = [ Coeff1(1) - Coeff1(2)*y(i,:) -
> > ) - Coeff1
> > (3)*y(i,:).*y(i+1,:)./(1 + Gain_comp*y(i+1,:));...
> > delta_z*y(i,:).*y(i+1,:) - (1 +
> > + Gain_comp*y(i+1,:)).*(y
> > (i+1,:) - P_in(tmp)); ...
> > Linw_en/2*delta_z*y(i,:) + (1 +
> > + Gain_comp*y(i+1,:)).*(y
> > (i+2,:) - Phi_in(tmp))];
> >
> > for i = 7:3:3*(N-1)-2
> > out(i-3:i-1,:) = [ Coeff1(1) - Coeff1(2)*y(i,:) -
> > ) - Coeff1
> > (3)*y(i,:).*y(i+1,:)./(1 + Gain_comp*y(i+1,:));...
> > delta_z*y(i,:).*y(i+1,:) - (1 +
> > 1 + Gain_comp*y(i+1,:)).*(y
> > (i+1,:) - y(i-2,:)); ...
> > Linw_en/2*delta_z*y(i,:) + (1 +
> > 1 + Gain_comp*y(i+1,:)).*(y
> > (i+2,:) - y(i-1,:))];
> > end
> > end
> >
> > The problem now is that output of the nested
> > nested function
> > (see varible 'out') is vector of 3(N-2)-by-1,while
> > the size
> > of y is 3*(N-1)-by-1. Thus it can not be performed
> > by
> > Matlab.
> >
> > Do I need to set the first 3 rows of 'out' to
> > ut' to be
> > zeros? But it seems not correct.
> >
> > Please help me have a look.
> > Many thanks.
> >
>
> It's just a question of enumeration of the variables.
> The system
>
> du1(i)/dt = c1 - c2*u1(i)-c3*u1(i)*u2(i)/(1+g*u2(i))
> 0 = (u2(i)-u2(i-1))/deltaz -u1(i)*u2(i)/(1+g*u2(i))
> 0 = (u3(i)-u3(i-1))/deltaz + l/2*u1(i)/(1+g*u2(i))
> (2 <= i <= n)
>
> has 3*(n-1) equations for the 3*(n-1) unknowns
> u1(2),...,u1(n),u2(2),...,u2(n),u3(2),...,u3(n).
>
> You have to shift the enumeration for one index
> (2 becomes 1, 3 becomes 2 etc.)
>
> Another way is to add the equations
> du1(1)/dt = c1 - c2*u1(1)-c3*u1(1)*u2(1)/(1+g*u2(1))
> 0 = u2(1) - P_in(tmp)
> 0 = u3(1) - Phi_in(tmp)
> and leave the enumeration as it is.
>
> Best wishes
> Torsten.

Hi, Torsten, thank you so much for all your patient and
clear explanations. Now the problem becomes "Need a better
guess y0 for consistent initial conditions.". should be
much smoother after I check the Matlab help manual later.

Best wishes and a good weekend!!

"Jinlong Wei" <jinlongwei@gmail.com> wrote in message
<g72a53$ktj$1@fred.mathworks.com>...
> Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote in
> message
>
<15852838.1217677313101.JavaMail.jakarta@nitrogen.mathforum.
> org>...
> > > Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote
> > > in
> > > message
> > > >
> > > > I did not check everything in detail, but we
> > > overlooked
> > > > the fact that our discretization gives no equation
> > > > to determine u3(z2).
> > > >
> > > > So maybe for a first test, a first-order
> > > approximation
> > > > for the spatial derivatives is sufficient:
> > > >
> > > > du1(i)/dt = c1 -
> > > c2*u1(i)-c3*u1(i)*u2(i)/(1+g*u2(i))
> > > > 0 = (u2(i)-u2(i-1))/deltaz -
> > > u1(i)*u2(i)/(1+g*u2(i))
> > > > 0 = (u3(i)-u3(i-1))/deltaz + l/2*u1(i)/(1+g*u2(i))
> > > > (2 <= i <= n)
> > > >
> > > > Best wishes
> > > > Torsten.
> > >
> > > Hi, Torsten,
> > > the new nested function is listed as following:
> > >
> > > function out = f(t, y)
> > > tmp = floor(t*fs)+1;
> > > %first-order approximation
> > > i = 4;
> > > out(i-3:i-1,:) = [ Coeff1(1) - Coeff1(2)*y(i,:) -
> > > ) - Coeff1
> > > (3)*y(i,:).*y(i+1,:)./(1 + Gain_comp*y(i+1,:));...
> > > delta_z*y(i,:).*y(i+1,:) - (1 +
> > > + Gain_comp*y(i+1,:)).*(y
> > > (i+1,:) - P_in(tmp)); ...
> > > Linw_en/2*delta_z*y(i,:) + (1 +
> > > + Gain_comp*y(i+1,:)).*(y
> > > (i+2,:) - Phi_in(tmp))];
> > >
> > > for i = 7:3:3*(N-1)-2
> > > out(i-3:i-1,:) = [ Coeff1(1) - Coeff1(2)*y(i,:) -
> > > ) - Coeff1
> > > (3)*y(i,:).*y(i+1,:)./(1 + Gain_comp*y(i+1,:));...
> > > delta_z*y(i,:).*y(i+1,:) - (1 +
> > > 1 + Gain_comp*y(i+1,:)).*(y
> > > (i+1,:) - y(i-2,:)); ...
> > > Linw_en/2*delta_z*y(i,:) + (1 +
> > > 1 + Gain_comp*y(i+1,:)).*(y
> > > (i+2,:) - y(i-1,:))];
> > > end
> > > end
> > >
> > > The problem now is that output of the nested
> > > nested function
> > > (see varible 'out') is vector of 3(N-2)-by-1,while
> > > the size
> > > of y is 3*(N-1)-by-1. Thus it can not be performed
> > > by
> > > Matlab.
> > >
> > > Do I need to set the first 3 rows of 'out' to
> > > ut' to be
> > > zeros? But it seems not correct.
> > >
> > > Please help me have a look.
> > > Many thanks.
> > >
> >
> > It's just a question of enumeration of the variables.
> > The system
> >
> > du1(i)/dt = c1 - c2*u1(i)-c3*u1(i)*u2(i)/(1+g*u2(i))
> > 0 = (u2(i)-u2(i-1))/deltaz -u1(i)*u2(i)/(1+g*u2(i))
> > 0 = (u3(i)-u3(i-1))/deltaz + l/2*u1(i)/(1+g*u2(i))
> > (2 <= i <= n)
> >
> > has 3*(n-1) equations for the 3*(n-1) unknowns
> > u1(2),...,u1(n),u2(2),...,u2(n),u3(2),...,u3(n).
> >
> > You have to shift the enumeration for one index
> > (2 becomes 1, 3 becomes 2 etc.)
> >
> > Another way is to add the equations
> > du1(1)/dt = c1 - c2*u1(1)-c3*u1(1)*u2(1)/(1+g*u2(1))
> > 0 = u2(1) - P_in(tmp)
> > 0 = u3(1) - Phi_in(tmp)
> > and leave the enumeration as it is.
> >
> > Best wishes
> > Torsten.
>
> Hi, Torsten, thank you so much for all your patient and
> clear explanations. Now the problem becomes "Need a
better
> guess y0 for consistent initial conditions.". should be
> much smoother after I check the Matlab help manual later.
>
> Best wishes and a good weekend!!

Hi, Torsten, I decided to solve u1 and u2 in a first step
because the consistent initial values can not be satisfied
if u3 is also included. The new problem for solving u1 and
u2 is the step size and tolerance as cited below:

Warning: Failure at t=2.328082e-012. Unable to meet
integration tolerances without reducing the step size below
the smallest value allowed (6.462349e-027) at time t.

I even reduced the AbsTol to 1E-3, still the same problem.
where is the smallest value allowed (6.462349e-027)
specified? can it be modified or not ? or any other
possible solutions?

Many thanks.
Jinlong

>Hi, Torsten, I decided to solve u1 and u2 in a first >step
>because the consistent initial values can not be >satisfied
>if u3 is also included. The new problem for solving u1 >and
>u2 is the step size and tolerance as cited below:
>
>Warning: Failure at t=2.328082e-012. Unable to meet
>integration tolerances without reducing the step size >below
>the smallest value allowed (6.462349e-027) at time t.
>
>I even reduced the AbsTol to 1E-3, still the same >problem.
>where is the smallest value allowed (6.462349e-027)
>specified? can it be modified or not ? or any other
>possible solutions?
>
>Many thanks.
>Jinlong

To see whether you programmed your problem correctly,
I'd suggest to first solve the following simpler
problem:

Du1/Dt = m - u1 - u1*u2/(1+r*u2)
0 = Du2/Dx - 1
0 = Du3/Dx + 1

Do you succeed in finding a solution ?

If consistent initial values are the problem, then
maybe you have the possibility to solve the nonlinear
system

0 = Du2/Dx - u1*u2/((1+r*u2))
0 = Du3/Dx + a*u1*u2/((1+r*u2))

at time t=0 for u2 and u3 by inserting the initial
values for u1 in the grid points.
The obtained values could then be used to call the
DAE solver with consistent values.

Best wishes
Torsten.

Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote in
message
<30145878.1217833527479.JavaMail.jakarta@nitrogen.mathforum.
org>...
> >Hi, Torsten, I decided to solve u1 and u2 in a first
>step
> >because the consistent initial values can not be
>satisfied
> >if u3 is also included. The new problem for solving u1
>and
> >u2 is the step size and tolerance as cited below:
> >
> >Warning: Failure at t=2.328082e-012. Unable to meet
> >integration tolerances without reducing the step size
>below
> >the smallest value allowed (6.462349e-027) at time t.
> >
> >I even reduced the AbsTol to 1E-3, still the same
>problem.
> >where is the smallest value allowed (6.462349e-027)
> >specified? can it be modified or not ? or any other
> >possible solutions?
> >
> >Many thanks.
> >Jinlong
>
> To see whether you programmed your problem correctly,
> I'd suggest to first solve the following simpler
> problem:
>
> Du1/Dt = m - u1 - u1*u2/(1+r*u2)
> 0 = Du2/Dx - 1
> 0 = Du3/Dx + 1
>
> Do you succeed in finding a solution ?

Hi, Torsten. I agree with you to try the simpler problem
first. In order to avoid the possible error from sparse
Jacobian matrix, I removed the 'JPattern' from options. The
result also show the same problem. that is,
Warning: Failure at t=4.650000e-012. Unable to meet
integration tolerances without reducing the step size below
the smallest value allowed (1.292470e-026) at time t.

I paste the key code here in case you want to have a look

   M = zeros(3*(N-1),3*(N-1));
   for ii = 1:3:3*(N-1)
       M(ii,ii) = 1;
   end
   val_AbsTol = 1e-10*ones(1,3*(N-1));
   tspan = 0: 1/fs: (length(Opt_in) - 1)/fs;
   y0 = repmat([G_zero; 0; 0],N-1,1);
   options = odeset('Mass',M,'absTol',
val_AbsTol,'Vectorized','on');
   
   [t,y] = ode15s(@f,tspan,y0,options);

    function out = f(t, y)
      tmp = floor(t*fs)+1;
      i = 1;
   out(i:i+2,:) = [ - Coeff1(3)*G_zero.*y(i+1,:)./(1 +
Gain_comp*y(i+1,:));...
                      y(i+1,:) - P_in(tmp); ...
                      y(i+2,:) - Phi_in(tmp)];
      for i = 4:3:3*(N-1)-2
   out(i:i+2,:) = [ Coeff1(1) - Coeff1(2)*y(i,:) - Coeff1(3)
*y(i,:).*y(i+1,:)./(1 + Gain_comp*y(i+1,:));...
                      delta_z - (y(i+1,:) - y(i-2,:)); ...
                      delta_z + (y(i+2,:) - y(i-1,:))]; %
first-order approximation
   end

Now,the simpler one even does not work. what should the
problem be ?

Many thanks.

>
> If consistent initial values are the problem, then
> maybe you have the possibility to solve the nonlinear
> system
>
> 0 = Du2/Dx - u1*u2/((1+r*u2))
> 0 = Du3/Dx + a*u1*u2/((1+r*u2))
>
> at time t=0 for u2 and u3 by inserting the initial
> values for u1 in the grid points.
> The obtained values could then be used to call the
> DAE solver with consistent values.
>
> Best wishes
> Torsten.

> Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote
> in
> message
> <30145878.1217833527479.JavaMail.jakarta@nitrogen.math
> forum.
> org>...
> Hi, Torsten. I agree with you to try the simpler
> problem
> first. In order to avoid the possible error from
> sparse
> Jacobian matrix, I removed the 'JPattern' from
> options. The
> result also show the same problem. that is,
> Warning: Failure at t=4.650000e-012. Unable to meet
> integration tolerances without reducing the step size
> below
> the smallest value allowed (1.292470e-026) at time t.
>
> I paste the key code here in case you want to have a
> look
>
> M = zeros(3*(N-1),3*(N-1));
> for ii = 1:3:3*(N-1)
> M(ii,ii) = 1;
> end
> val_AbsTol = 1e-10*ones(1,3*(N-1));
> tspan = 0: 1/fs: (length(Opt_in) - 1)/fs;
> y0 = repmat([G_zero; 0; 0],N-1,1);
> options = odeset('Mass',M,'absTol',
> val_AbsTol,'Vectorized','on');
>
> [t,y] = ode15s(@f,tspan,y0,options);
>
> function out = f(t, y)
> tmp = floor(t*fs)+1;
> i = 1;
> out(i:i+2,:) = [ - Coeff1(3)*G_zero.*y(i+1,:)./(1
> /(1 +
> Gain_comp*y(i+1,:));...
> y(i+1,:) - P_in(tmp); ...
> y(i+2,:) - Phi_in(tmp)];
> for i = 4:3:3*(N-1)-2
> out(i:i+2,:) = [ Coeff1(1) - Coeff1(2)*y(i,:) -
> ) - Coeff1(3)
> *y(i,:).*y(i+1,:)./(1 + Gain_comp*y(i+1,:));...
> delta_z - (y(i+1,:) -
> delta_z - (y(i+1,:) - y(i-2,:)); ...
> delta_z + (y(i+2,:) -
> delta_z + (y(i+2,:) - y(i-1,:))]; %
> first-order approximation
> end
>
> Now,the simpler one even does not work. what should
> the
> problem be ?
>
> Many thanks.
>
>

I suspect that your integration is not successful
because you don't supply consistent initial values
for u2 and u3.
For the test problem, use the following consistent
values as initial values for the variables:

u1(i) = u1(i,t=0) (I think these values are in the
G_zero array ?)
u2(i) = P_in + (i-1)*deltaz (i=1,...,N-1)
u3(i) = Phi_in - (i-1)*deltaz (i=1,...,N-1)

Is y0 your vector of initial values for the variables ?
Then beware that they are arranged in y0 in the same
order as in the vector y used in the function out
(I think this is not the case in your matlab code):

y0(3*i-2) contains u1(i,t=0)
y0(3*i-1) contains u2(i,t=0)
y0(3*i) contains u3(i,t=0)

With these modification, the solver should start
the integration.

If not, you should begin with a very simple problem,
e.g.
du1/dt = u1
0 = u2 - u1

Best wishes
Torsten.

> > Torsten Hennig <Torsten.Hennig@umsicht.fhg.de>
> wrote
> > in
> > message
> >
> <30145878.1217833527479.JavaMail.jakarta@nitrogen.math
>
> > forum.
> > org>...
> I suspect that your integration is not successful
> because you don't supply consistent initial values
> for u2 and u3.
> For the test problem, use the following consistent
> values as initial values for the variables:
>
> u1(i) = u1(i,t=0) (I think these values are in the
> G_zero array ?)
> u2(i) = P_in + (i-1)*deltaz (i=1,...,N-1)
> u3(i) = Phi_in - (i-1)*deltaz (i=1,...,N-1)
>
> Is y0 your vector of initial values for the variables
> ?
> Then beware that they are arranged in y0 in the same
> order as in the vector y used in the function out
> (I think this is not the case in your matlab code):
>
> y0(3*i-2) contains u1(i,t=0)
> y0(3*i-1) contains u2(i,t=0)
> y0(3*i) contains u3(i,t=0)

I see now that the order of the variables is correct,
but you should make an attempt to supply the indicated
consistent initial values.

The following won't be the reason why the solver
doesn't start the interation, but why do you supply a different ODE for u1(1) ?
Why not setting
du1(1)/dt= Coeff1(1) - Coeff1(2)*u1(1) -
 - Coeff1(3)*u1(1)*u2(1)/(1 + Gain_comp*u2(1)) ?
as for the other grid points ?

>
> With these modification, the solver should start
> the integration.
>
> If not, you should begin with a very simple problem,
> e.g.
> du1/dt = u1
> 0 = u2 - u1
>
> Best wishes
> Torsten.

Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote in
message
<25427037.1217918854575.JavaMail.jakarta@nitrogen.mathforum.
org>...
> > Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote
> > in
> > message
> > <30145878.1217833527479.JavaMail.jakarta@nitrogen.math
> > forum.
> > org>...
> > Hi, Torsten. I agree with you to try the simpler
> > problem
> > first. In order to avoid the possible error from
> > sparse
> > Jacobian matrix, I removed the 'JPattern' from
> > options. The
> > result also show the same problem. that is,
> > Warning: Failure at t=4.650000e-012. Unable to meet
> > integration tolerances without reducing the step size
> > below
> > the smallest value allowed (1.292470e-026) at time t.
> >
> > I paste the key code here in case you want to have a
> > look
> >
> > M = zeros(3*(N-1),3*(N-1));
> > for ii = 1:3:3*(N-1)
> > M(ii,ii) = 1;
> > end
> > val_AbsTol = 1e-10*ones(1,3*(N-1));
> > tspan = 0: 1/fs: (length(Opt_in) - 1)/fs;
> > y0 = repmat([G_zero; 0; 0],N-1,1);
> > options = odeset('Mass',M,'absTol',
> > val_AbsTol,'Vectorized','on');
> >
> > [t,y] = ode15s(@f,tspan,y0,options);
> >
> > function out = f(t, y)
> > tmp = floor(t*fs)+1;
> > i = 1;
> > out(i:i+2,:) = [ - Coeff1(3)*G_zero.*y(i+1,:)./(1
> > /(1 +
> > Gain_comp*y(i+1,:));...
> > y(i+1,:) - P_in(tmp); ...
> > y(i+2,:) - Phi_in(tmp)];
> > for i = 4:3:3*(N-1)-2
> > out(i:i+2,:) = [ Coeff1(1) - Coeff1(2)*y(i,:) -
> > ) - Coeff1(3)
> > *y(i,:).*y(i+1,:)./(1 + Gain_comp*y(i+1,:));...
> > delta_z - (y(i+1,:) -
> > delta_z - (y(i+1,:) - y(i-2,:)); ...
> > delta_z + (y(i+2,:) -
> > delta_z + (y(i+2,:) - y(i-1,:))]; %
> > first-order approximation
> > end
> >
> > Now,the simpler one even does not work. what should
> > the
> > problem be ?
> >
> > Many thanks.
> >
> >
>
> I suspect that your integration is not successful
> because you don't supply consistent initial values
> for u2 and u3.
> For the test problem, use the following consistent
> values as initial values for the variables:
>
> u1(i) = u1(i,t=0) (I think these values are in the
> G_zero array ?)
> u2(i) = P_in + (i-1)*deltaz (i=1,...,N-1)
> u3(i) = Phi_in - (i-1)*deltaz (i=1,...,N-1)
>
> Is y0 your vector of initial values for the variables ?
> Then beware that they are arranged in y0 in the same
> order as in the vector y used in the function out
> (I think this is not the case in your matlab code):
>
> y0(3*i-2) contains u1(i,t=0)
> y0(3*i-1) contains u2(i,t=0)
> y0(3*i) contains u3(i,t=0)
>
> With these modification, the solver should start
> the integration.
>
> If not, you should begin with a very simple problem,
> e.g.
> du1/dt = u1
> 0 = u2 - u1

Hi, Torsten, I was absent for a period of time and now
continue the work. hope you will still give me suggestions.
I tried a very simple one:
du1/dt = a*u1
     0 = u2 - u1
     0 = u3 + u1
It works fine.

But when I try the following one as you suggested:
Du1/Dt = m - u1 - u1*u2/(1+r*u2)
0 = Du2/Dx - 1
0 = Du3/Dx + 1
and I also set initial value of u2 and u3 as
u2(i) = p_in + (i-1)*delta_z
u3(i) = Phi_in + (i-1)*delta_z
which I think are consistent.

But the integration still failed. It still stops at a point
and complains about the step size can not meet integration
tolerance even the step size is the smallest allowed.

Is there a convinient and fast way to set and check the
initial values to be consistent?

Thanks and best wishes.

>
> Best wishes
> Torsten.

> Hi, Torsten, I was absent for a period of time and
> now
> continue the work. hope you will still give me
> suggestions.
> I tried a very simple one:
> du1/dt = a*u1
> 0 = u2 - u1
> 0 = u3 + u1
> It works fine.
>
> But when I try the following one as you suggested:
> Du1/Dt = m - u1 - u1*u2/(1+r*u2)
> 0 = Du2/Dx - 1
> 0 = Du3/Dx + 1
> and I also set initial value of u2 and u3 as
> u2(i) = p_in + (i-1)*delta_z
> u3(i) = Phi_in + (i-1)*delta_z

u3(i) = Phi_in - (i-1)*delta_z !!

> which I think are consistent.
>
> But the integration still failed. It still stops at a
> point
> and complains about the step size can not meet
> integration
> tolerance even the step size is the smallest allowed.
>
> Is there a convinient and fast way to set and check
> the
> initial values to be consistent?
>

Yes.
If it is difficult to compute the algebraic variables
at time t = 0 analytically, you can do the following:

Fix the differential variables to their initial values, take a guess for the algebraic variables and solve the
part of the DAE-system which contains the algebraic
equations with fsolve.
This should give you a vector of starting values for
the full DAE-system which is consistent.

In the case above, at t=0, fix u1(i), take a guess for
u2(i) and u3(i) and solve the 2*(N-1)
algebraic equations for u2(i) and u3(i) with fsolve.

> Thanks and best wishes.
>
> >

Best wishes
Torsten.

Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote in
message
<9732519.1219302095610.JavaMail.jakarta@nitrogen.mathforum.o
rg>...
> > Hi, Torsten, I was absent for a period of time and
> > now
> > continue the work. hope you will still give me
> > suggestions.
> > I tried a very simple one:
> > du1/dt = a*u1
> > 0 = u2 - u1
> > 0 = u3 + u1
> > It works fine.
> >
> > But when I try the following one as you suggested:
> > Du1/Dt = m - u1 - u1*u2/(1+r*u2)
> > 0 = Du2/Dx - 1
> > 0 = Du3/Dx + 1
> > and I also set initial value of u2 and u3 as
> > u2(i) = p_in + (i-1)*delta_z
> > u3(i) = Phi_in + (i-1)*delta_z
>
> u3(i) = Phi_in - (i-1)*delta_z !!
>
> > which I think are consistent.
> >
> > But the integration still failed. It still stops at a
> > point
> > and complains about the step size can not meet
> > integration
> > tolerance even the step size is the smallest allowed.
> >
> > Is there a convinient and fast way to set and check
> > the
> > initial values to be consistent?
> >
>
> Yes.
> If it is difficult to compute the algebraic variables
> at time t = 0 analytically, you can do the following:
>
> Fix the differential variables to their initial values,
take a guess for the algebraic variables and solve the
> part of the DAE-system which contains the algebraic
> equations with fsolve.
> This should give you a vector of starting values for
> the full DAE-system which is consistent.
>
> In the case above, at t=0, fix u1(i), take a guess for
> u2(i) and u3(i) and solve the 2*(N-1)
> algebraic equations for u2(i) and u3(i) with fsolve.

Hi, Torsten. I did exactly as you suggested. and it works
for
 Du1/Dt = m - u1 - u1*u2/(1+r*u2)
 0 = Du2/Dx - 1
 0 = Du3/Dx + 1
I also checked that the exitflag = 1 which means fsolve
succeds. and the integration of ode is also successful.

The problem is, when I do the same thing to
 Du1/Dt = m - u1 - u1*u2/(1+r*u2)
      0 = Du2/Dx - u1*u2/(1+r*u2)
      0 = Du3/Dx + c*u1/(1+r*u2)
that is, fix u1(i) to thier initial values and solve the 2*
(N-1) algebraic equations for u2(i) and u3(i) with fsolve.

   I checked that fsolve works because the exitflag is 1.
which means the initial values are consistent now.

   But the integration still has the same problem:
stops at a point and complains about the step size can not
meet integration tolerance even the step size is the
smallest allowed.

   Now it seems very close to success. what do you think
should be the problem in this step ?

   I'm looking forward to your help.
   Many thanks and wishes.

> Hi, Torsten. I did exactly as you suggested. and it
> works
> for
> Du1/Dt = m - u1 - u1*u2/(1+r*u2)
> 0 = Du2/Dx - 1
> 0 = Du3/Dx + 1
> I also checked that the exitflag = 1 which means
> fsolve
> succeds. and the integration of ode is also
> successful.
>
> The problem is, when I do the same thing to
> Du1/Dt = m - u1 - u1*u2/(1+r*u2)
> 0 = Du2/Dx - u1*u2/(1+r*u2)
> 0 = Du3/Dx + c*u1/(1+r*u2)
> that is, fix u1(i) to thier initial values and solve
> the 2*
> (N-1) algebraic equations for u2(i) and u3(i) with
> fsolve.
>
> I checked that fsolve works because the exitflag
> lag is 1.
> which means the initial values are consistent now.
>
> But the integration still has the same problem:
> stops at a point and complains about the step size
> can not
> meet integration tolerance even the step size is the
>
> smallest allowed.
>
> Now it seems very close to success. what do you
> you think
> should be the problem in this step ?
>
> I'm looking forward to your help.
> Many thanks and wishes.
>

If you use the approximation
du2(i)/dx ~ (u2(i)-u2(i-1))/deltaz,
the equation
0 = Du2/Dx - u1*u2/(1+r*u2)
results in
(u2(i)-u2(i-1))/deltaz - u1(i)*u2(i)/(1+r*u2(i)).
Solving for u2(i) gives
u2(i) = -(1-deltaz*u1(i)-r*u2(i-1))/(2*r) + or -
sqrt(((1-deltaz*u1(i)-r*u2(i-1))/(2*r))^2 + u2(i-1)/r).
You see the problem:
((1-deltaz*u1(i)-r*u2(i-1))/(2*r))^2 + u2(i-1)/r
may become negative during integration so that the
sqrt can not be taken, or you may get problems
because there are two possible values for u2
both satisfying your second equation.

Best wishes
Torsten.

Torsten Hennig <Torsten.Hennig@umsicht.fhg.de> wrote in
message
<29168110.1219647213478.JavaMail.jakarta@nitrogen.mathforum.
org>...
> > Hi, Torsten. I did exactly as you suggested. and it
> > works
> > for
> > Du1/Dt = m - u1 - u1*u2/(1+r*u2)
> > 0 = Du2/Dx - 1
> > 0 = Du3/Dx + 1
> > I also checked that the exitflag = 1 which means
> > fsolve
> > succeds. and the integration of ode is also
> > successful.
> >
> > The problem is, when I do the same thing to
> > Du1/Dt = m - u1 - u1*u2/(1+r*u2)
> > 0 = Du2/Dx - u1*u2/(1+r*u2)
> > 0 = Du3/Dx + c*u1/(1+r*u2)
> > that is, fix u1(i) to thier initial values and solve
> > the 2*
> > (N-1) algebraic equations for u2(i) and u3(i) with
> > fsolve.
> >
> > I checked that fsolve works because the exitflag
> > lag is 1.
> > which means the initial values are consistent now.
> >
> > But the integration still has the same problem:
> > stops at a point and complains about the step size
> > can not
> > meet integration tolerance even the step size is the
> >
> > smallest allowed.
> >
> > Now it seems very close to success. what do you
> > you think
> > should be the problem in this step ?
> >
> > I'm looking forward to your help.
> > Many thanks and wishes.
> >
>
> If you use the approximation
> du2(i)/dx ~ (u2(i)-u2(i-1))/deltaz,
> the equation
> 0 = Du2/Dx - u1*u2/(1+r*u2)
> results in
> (u2(i)-u2(i-1))/deltaz - u1(i)*u2(i)/(1+r*u2(i)).
> Solving for u2(i) gives
> u2(i) = -(1-deltaz*u1(i)-r*u2(i-1))/(2*r) + or -
> sqrt(((1-deltaz*u1(i)-r*u2(i-1))/(2*r))^2 + u2(i-1)/r).
> You see the problem:
> ((1-deltaz*u1(i)-r*u2(i-1))/(2*r))^2 + u2(i-1)/r
> may become negative during integration so that the
> sqrt can not be taken, or you may get problems
> because there are two possible values for u2
> both satisfying your second equation.

  Actually, what you said may not be the problem. Because I
checked the values of u2 for those successfully integrated
points before the stop point, I found all of them are
positve values and u2 increases as time steps forward.

  I also checked the calculation of consistent values using
fsolve. I found that though the exitflag = 1, there is
still a diagnositc message, that is

 Optimization terminated: first-order optimality is less
than options.TolFun.

  I think the initial values for u1,u2 and u3 are still not
consistent. Because I only calculated u2(i) and u3(i) by
using fsolve, while u1(i) are fixed. Physically the fixed u
(i) are correct, but in this situation, together with
initial values of u2(i) and u3(i), they are not consistent
any more.

  I changed the F in the nested function of fsolve. that
is,
 0 = m - u1 - u1*u2/(1+r*u2)
 0 = Du2/Dx - u1*u2/(1+r*u2)
 0 = Du3/Dx + c*u1/(1+r*u2)
  to solve the initial values together. But the calculation
seems never ends.

  Thanks and best wishes.
  Jinlong

> > If you use the approximation
> > du2(i)/dx ~ (u2(i)-u2(i-1))/deltaz,
> > the equation
> > 0 = Du2/Dx - u1*u2/(1+r*u2)
> > results in
> > (u2(i)-u2(i-1))/deltaz - u1(i)*u2(i)/(1+r*u2(i)).
> > Solving for u2(i) gives
> > u2(i) = -(1-deltaz*u1(i)-r*u2(i-1))/(2*r) + or -
> > sqrt(((1-deltaz*u1(i)-r*u2(i-1))/(2*r))^2 +
> u2(i-1)/r).
> > You see the problem:
> > ((1-deltaz*u1(i)-r*u2(i-1))/(2*r))^2 + u2(i-1)/r
> > may become negative during integration so that the
> > sqrt can not be taken, or you may get problems
> > because there are two possible values for u2
> > both satisfying your second equation.
>
> Actually, what you said may not be the problem.
> m. Because I
> checked the values of u2 for those successfully
> integrated
> points before the stop point, I found all of them are
>
> positve values and u2 increases as time steps
> forward.
>
> I also checked the calculation of consistent values
> es using
> fsolve. I found that though the exitflag = 1, there
> is
> still a diagnositc message, that is
>
> Optimization terminated: first-order optimality is
> s less
> than options.TolFun.
>
> I think the initial values for u1,u2 and u3 are
> re still not
> consistent. Because I only calculated u2(i) and u3(i)
> by
> using fsolve, while u1(i) are fixed. Physically the
> fixed u
> (i) are correct, but in this situation, together with
>
> initial values of u2(i) and u3(i), they are not
> consistent
> any more.
>
> I changed the F in the nested function of fsolve.
> e. that
> is,
> 0 = m - u1 - u1*u2/(1+r*u2)
> 0 = Du2/Dx - u1*u2/(1+r*u2)
> 0 = Du3/Dx + c*u1/(1+r*u2)

You must not include the first equation for the
calculation of consistent values by fsolve because
this is the differential one.
Only take the second and third equations with
u1(i) = u1(i,t=0) fixed (initial values
for u1 at time t=0).

> to solve the initial values together. But the
> he calculation
> seems never ends.
>
> Thanks and best wishes.
> Jinlong
>
>

If the above correction also leads to nowhere, you should
use a normal ODE - integrator for u1 and calculate the
values for u2 and u3 explicitly at each call for the
right hand sides.

u2(1) and u3(1) are always given as boundary
conditions. Now with the above recursion

u2(i) = -(1-deltaz*u1(i)-r*u2(i-1))/(2*r) + or -
 sqrt(((1-deltaz*u1(i)-r*u2(i-1))/(2*r))^2 +
 u2(i-1)/r),

you can calculate u2(i) _explicitly_ at each time t
given the value for u1(i) from the integrator.
(Of course, you must know from the physics behind the
equations whether to add or subtract the sqrt).

Now you have u1(i) and u2(i)for all i.

With this information, u3(i) can also be calculated
recursively:

(u3(i)-u3(i-1))/deltaz + c*u1(i)/(1+r*u2(i)) ->
u3(i) = u3(i-1) - deltaz*c*u1(i)/(1+r*u2(i))

This method should work without any problems.

Best wishes
Torsten.