function writing in differential equations

[genehacker]

Hi,

     I am using dassl to solve algebraic differential equations and want to write a variable as a function of time according to the dassl condition 0 = f(x,xdot,t)

If I tried writing x(1) as a function only on time exponentially, like this, x(1) = exp(-10*t), it gives an error. I would like to know the correct way of writing a variable dependent only on time, if the function that solves the ADE is function res = f(x,xdot,t);
Please help.

ps: I made sure that the initial conditions for x and xdot are correct, I dont think the error is due to that.

cheers
k3

This is the error generated
----------------------------
 ***MESSAGE FROM ROUTINE DDASSL IN LIBRARY SLATEC.
 ***POTENTIALLY RECOVERABLE ERROR, PROG ABORTED, TRACEBACK REQUESTED
 *  AT T =    0.000000E+00  AND STEPSIZE H =    3.215509E-16  THE ERROR TEST
 *  FAILED REPEATEDLY OR WITH ABS(H)=HMIN
 *  ERROR NUMBER = -6
 *
 ***END OF MESSAGE

 ***JOB ABORT DUE TO UNRECOVERED ERROR.
0          ERROR MESSAGE SUMMARY
 LIBRARY    SUBROUTINE MESSAGE START             NERR     LEVEL     COUNT
 SLATEC     DDASSL     AT CURRENT T =    1.        -1         1         2
 SLATEC     DDASSL     AT T =    0.000000E+        -8         1        12
 SLATEC     DDASSL     AT T =    0.000000E+        -6         1        25
 SLATEC     DDASSL     AT T =    4.710277E-        -6         1         2
 SLATEC     DDASSL     AT T =    1.538370E-        -6         1         2
 SLATEC     DDASSL     AT T =    1.655685E-        -6         1         3
 SLATEC     DDASSL     AT T =    0.000000E+        -7         1         4
 SLATEC     DDASSL     AT T =    2.000000E-        -3         1         2
 SLATEC     DDASSL     AT T =    3.762919E-        -6         1         7
 SLATEC     DDASSL     AT T =    1.505168E-        -6         1         2
0OTHER ERRORS NOT INDIVIDUALLY TABULATED =         12

[Olaf Till]

On Mon, Oct 13, 2008 at 04:06:26PM -0700, genehacker wrote:

>
> Hi,
>
>      I am using dassl to solve algebraic differential equations and want to
> write a variable as a function of time according to the dassl condition 0 =
> f(x,xdot,t)
>
> If I tried writing x(1) as a function only on time exponentially, like this,
> x(1) = exp(-10*t), it gives an error. I would like to know the correct way
> of writing a variable dependent only on time, if the function that solves
> the ADE is function res = f(x,xdot,t);
> Please help.

A variable which only depends on time does not belong into the vector
'x'. Cancel the equation for 'x(1)' from your function 'f'. If some
other equation in 'f' needed 'x(1)', substitute 'exp(-10*t)'
instead. If some other equation in 'f' needed 'xdot(1)', substitute
the derivative of 'x(1)' (-10*exp(-10*t)) instead.

Any variable which does not depend only on time, but its derivative
plays no role and the variable is not needed by another equation in
'f', should be calculated from the values of 'x', 'xdot', and 't'
after the integration is ready.

You should consider this and check if you really need a DAE solver,
not an ODE solver (which a former posting of yours suggested).

HTH, Olaf
_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[Carlo de Falco-2]

On 14/ott/08, at 09:11, Olaf Till wrote:

> On Mon, Oct 13, 2008 at 04:06:26PM -0700, genehacker wrote:
>>
>> Hi,
>>
>>     I am using dassl to solve algebraic differential equations and  
>> want to
>> write a variable as a function of time according to the dassl  
>> condition 0 =
>> f(x,xdot,t)
>>
>> If I tried writing x(1) as a function only on time exponentially,  
>> like this,
>> x(1) = exp(-10*t), it gives an error. I would like to know the  
>> correct way
>> of writing a variable dependent only on time, if the function that  
>> solves
>> the ADE is function res = f(x,xdot,t);
>> Please help.
>
> A variable which only depends on time does not belong into the vector
> 'x'. Cancel the equation for 'x(1)' from your function 'f'. If some
> other equation in 'f' needed 'x(1)', substitute 'exp(-10*t)'
> instead. If some other equation in 'f' needed 'xdot(1)', substitute
> the derivative of 'x(1)' (-10*exp(-10*t)) instead.
>
> Any variable which does not depend only on time, but its derivative
> plays no role and the variable is not needed by another equation in
> 'f', should be calculated from the values of 'x', 'xdot', and 't'
> after the integration is ready.
>
> You should consider this and check if you really need a DAE solver,
> not an ODE solver (which a former posting of yours suggested).
>
> HTH, Olaf

genehacker,

If one component of the vector x is a known function of time there is  
no need to include it among the unknowns as Olaf noted.
I think there is rather a problem with mathematical notation than with  
Octave syntax here. Could you write in mathematical terms the problem  
you need to solve?

c.

_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[genehacker]

Thanx for the replies. So I directly put in the functional form in other equations as it has to be as you mentioned.
Now I dont understand whether I have to use DAE or ODE solver. I have been thinking of DAE because I have 'time-dependent variables' that are a function of other variables in the model,

So, If I have the following system of 6 equations:

xdot(1) = 0.5 * exp(-10*t) - x(4) * x(1);
xdot(2) = x(5) * 0.5 - x(6) * x(2);
xdot(3) = x(6) * x(2) - 0.5 * x(3);
x(4) = 0.5- ((0.5- 0.1)/(1 + (x(1)/5)));
x(5) = 0.5- ((0.5- 0.1)/(1 + (x(1)/5)));
x(6) = 0.1-((0.1- 0.5)/(1 + (x(1)/5)));

If x(4) = f(x(1)) is regarded as an algebraic equation, and x(4) is time-dependent and influences x(1), I would think there are two approaches to solve the problem and I have two questions in them.

1) Treat the system as DAE and use dassl to solve. But it requires to specify xdot0! Do we need to specify this at all? Isn't it over-solving the problem with more constraints? Is there a way not to specify x0 and still use dassl?

2) Only treat x(1),x(2) and x(3) as variables of the system and consider x(4),x(5) and x(6) as time-dependent parameters and pass these parameters to ODE at each time step. So it looks like a loop.

foreach timestep
{
      compute paramteres
      feed to ODE
}
In this case, how would I send these parameters to this function Lsode?

Correct me if Im wrong or doesnt make any sense.

Ps: I dont care about derivates of x(4),x(5) and x(6)

thanks
k3

[John W. Eaton]

On 14-Oct-2008, genehacker wrote:

|
| Thanx for the replies. So I directly put in the functional form in other
| equations as it has to be as you mentioned.
| Now I dont understand whether I have to use DAE or ODE solver. I have been
| thinking of DAE because I have 'time-dependent variables' that are a
| function of other variables in the model,
|
| So, If I have the following system of 6 equations:
|
| xdot(1) = 0.5 * exp(-10*t) - x(4) * x(1);
| xdot(2) = x(5) * 0.5 - x(6) * x(2);
| xdot(3) = x(6) * x(2) - 0.5 * x(3);
| x(4) = 0.5- ((0.5- 0.1)/(1 + (x(1)/5)));
| x(5) = 0.5- ((0.5- 0.1)/(1 + (x(1)/5)));
| x(6) = 0.1-((0.1- 0.5)/(1 + (x(1)/5)));
|
| If x(4) = f(x(1)) is regarded as an algebraic equation, and x(4) is
| time-dependent and influences x(1), I would think there are two approaches
| to solve the problem and I have two questions in them.
|
| 1) Treat the system as DAE and use dassl to solve. But it requires to
| specify xdot0! Do we need to specify this at all? Isn't it over-solving the
| problem with more constraints? Is there a way not to specify x0 and still
| use dassl?

To write your set of equations above in a form suitable for dassl (or
probably better, daspk, which is a later revision of dassl) you will
need to write them as:

  res(1) = xdot(1) - 0.5 * exp(-10*t) - x(4) * x(1);
  res(2) = xdot(2) - x(5) * 0.5 - x(6) * x(2);
  res(3) = xdot(3) - x(6) * x(2) - 0.5 * x(3);
  res(4) = x(4) - 0.5- ((0.5- 0.1)/(1 + (x(1)/5)));
  res(5) = x(5) - 0.5- ((0.5- 0.1)/(1 + (x(1)/5)));
  res(6) = x(6) - 0.1-((0.1- 0.5)/(1 + (x(1)/5)));

What is your initial condition?  You are correct that for best
results, it should be consistent, and dassl tends to have difficulty
if it is not.  If you use daspk, it has some options for attempting to
compute a consisten set of initial conditions (see "help
daspk_options" for details).  For example, I think you would need to
set

  daspk_options ("compute consistent initial condition", 1);
  daspk_options ("algebraic variables", [0, 0, 0, 1, 1, 1]);

before calling daspk.  With these settings, daspk can solve this
problem

  t = 0:0.1:10;
  x0 = zeros (6, 1);
  xdot0 = zeros (6, 1);

  [x, xdot, istate, msg] = daspk (@rhs, x0, xdot0, t);

and you can plot it with

  plot (t, x)

though your system of equations appears to be unstable (at least from
this initial condition).

| Ps: I dont care about derivates of x(4),x(5) and x(6)

You'll get them as part of the solution to the DAE system.

jwe
_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[Carlo de Falco-2]

On 14/ott/08, at 18:52, genehacker wrote:

> xdot(1) = 0.5 * exp(-10*t) - x(4) * x(1);
> xdot(2) = x(5) * 0.5 - x(6) * x(2);
> xdot(3) = x(6) * x(2) - 0.5 * x(3);
> x(4) = 0.5- ((0.5- 0.1)/(1 + (x(1)/5)));
> x(5) = 0.5- ((0.5- 0.1)/(1 + (x(1)/5)));
> x(6) = 0.1-((0.1- 0.5)/(1 + (x(1)/5)));

put in this form the system is an ODE as you can eliminate x(4:6) in  
terms of x(1).
You can solve the system this way:

x4 = @(x1) 0.5- ((0.5- 0.1)/(1 + (x1/5)));
x5 = @(x1) 0.5- ((0.5- 0.1)/(1 + (x1/5)));
x6 = @(x1) 0.1- ((0.1- 0.5)/(1 + (x1/5)));

xdot = @(x,t) [ 0.5 * exp(-10*t) - x4(x(1)) * x(1);
               x5(x(1)) * 0.5 - x6(x(1)) * x(2);
               x6(x(1)) * x(2) - 0.5 * x(3)];

t = linspace(0,10,100);
x = lsode (xdot, zeros(3,1), t);

plot(t,x)



c.

_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[Carlo de Falco-2]

On 14/ott/08, at 20:27, Carlo de Falco wrote:

>
> On 14/ott/08, at 18:52, genehacker wrote:
>
>> xdot(1) = 0.5 * exp(-10*t) - x(4) * x(1);
>> xdot(2) = x(5) * 0.5 - x(6) * x(2);
>> xdot(3) = x(6) * x(2) - 0.5 * x(3);
>> x(4) = 0.5- ((0.5- 0.1)/(1 + (x(1)/5)));
>> x(5) = 0.5- ((0.5- 0.1)/(1 + (x(1)/5)));
>> x(6) = 0.1-((0.1- 0.5)/(1 + (x(1)/5)));
>
> put in this form the system is an ODE as you can eliminate x(4:6) in  
> terms of x(1).
> You can solve the system this way:
>
> x4 = @(x1) 0.5- ((0.5- 0.1)/(1 + (x1/5)));
> x5 = @(x1) 0.5- ((0.5- 0.1)/(1 + (x1/5)));
> x6 = @(x1) 0.1- ((0.1- 0.5)/(1 + (x1/5)));
>
> xdot = @(x,t) [ 0.5 * exp(-10*t) - x4(x(1)) * x(1);
>              x5(x(1)) * 0.5 - x6(x(1)) * x(2);
>              x6(x(1)) * x(2) - 0.5 * x(3)];
>
> t = linspace(0,10,100);
> x = lsode (xdot, zeros(3,1), t);
>
> plot(t,x)
> c.

BTW if for some reason you still wish to use DASSL (or DASPK) for  
solving the system above, you can do:

f = @(x,xprime,t) ( xprime - xdot(x,t));
X = daspk (f, zeros(3,1), xdot(zeros(3,1),0), t);

which of course gives the same result as you can see by doing

plot(t,x,t,X,'kx')

c.



_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[boazilan]

Hello,

I am trying to install ARPACK for Octave.app on a Mac (G5 PPC Tiger).

I am following the steps in this thread:
https://www-old.cae.wisc.edu/pipermail/help-octave/2007-November.txt
I downloaded arpack-1.0.6.tar.gz from octave-forge and unpacked it to  
~/octave/arpack-1.0.6
I downloaded the ARPACK source from http://www.caam.rice.edu/software/ 
ARPACK/, compiled it, and moved libarpack.a to ~/octave/arpack-1.0.6/
In the subsequent make process, arpack-1.0.6/inst/svds.m calls for  
eigs.m, which I gather is supposed to be created using eigs.oct.

According to the thread above for arpack-1.0.2, eigs.oct is under  
powerpc-apple-darwin7.9.1-api-v29.
arpack-1.0.2 had the following directory structure:
 >  [~/octave/arpack-1.0.2] -bash-2.05b 517$ ls
 > doc                                     powerpc-apple-darwin7.9.1-
api-v29
 > packinfo                                svds.m

However, arpack-1.0.6 has a different directory structure:
% du -d1
48      ./doc
16      ./inst
3232    ./src

Apparently powerpc-apple-darwin7.9.1-api-v29 disappeared and with it  
eigs.oct.

Is this a problem in the octave-forge distribution package?

Boaz


_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[Carlo de Falco-2]

On 16/ott/08, at 02:47, Ilan Boaz wrote:

>
> Hello,
>
> I am trying to install ARPACK for Octave.app on a Mac (G5 PPC Tiger).
>
> I am following the steps in this thread:
> https://www-old.cae.wisc.edu/pipermail/help-octave/2007-November.txt
> I downloaded arpack-1.0.6.tar.gz from octave-forge and unpacked it to
> ~/octave/arpack-1.0.6
> I downloaded the ARPACK source from http://www.caam.rice.edu/software/
> ARPACK/, compiled it, and moved libarpack.a to ~/octave/arpack-1.0.6/
> In the subsequent make process, arpack-1.0.6/inst/svds.m calls for
> eigs.m, which I gather is supposed to be created using eigs.oct.
>
> According to the thread above for arpack-1.0.2, eigs.oct is under
> powerpc-apple-darwin7.9.1-api-v29.
> arpack-1.0.2 had the following directory structure:
>> [~/octave/arpack-1.0.2] -bash-2.05b 517$ ls
>> doc                                     powerpc-apple-darwin7.9.1-
> api-v29
>> packinfo                                svds.m
>
> However, arpack-1.0.6 has a different directory structure:
> % du -d1
> 48      ./doc
> 16      ./inst
> 3232    ./src
>
> Apparently powerpc-apple-darwin7.9.1-api-v29 disappeared and with it
> eigs.oct.
>
> Is this a problem in the octave-forge distribution package?
>
> Boaz
>

read the FAQ for the package manager at
http://octave.sourceforge.net/FAQ.html

or just type
 >> doc pkg
at the Octave prompt

summarizing: you don't just copy manually the files from a package to  
the install location, you need to install them using the package  
manager pkg.m

c.


_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[Carlo de Falco-2]

On 16/ott/08, at 17:29, Ilan Boaz wrote:

>>
>> summarizing: you don't just copy manually the files from a package  
>> to the install location, you need to install them using the package  
>> manager pkg.m
>
> That was the first thing I tried, but it breaks -- see below. It  
> seems to be looking for eigs.oct, which was why I unpacked the  
> package to look for it.

it tries to produce eigs.oct from the source code which is in the  
package and fails because it does not know where libarpack.a is as  
indicated by the error message:

> "/usr/libexec/gcc/powerpc-apple-darwin8/4.0.1/ld: can't locate file  
> for: -larpack
> collect2: ld returned 1 exit status"

  so no eigs.oct is created

> The thread
>>> https://www-old.cae.wisc.edu/pipermail/help-octave/2007-November.txt
> also suggests that eigs.oct should be inside this package.

eigs.oct is not in the package tarball it will be in the package  
installation directory after compilation but you won't get it there
until you manage to build it

> Any ideas?

Set the environment variable LDFLAGS to point to
LDFLAGS="-Lthe/directory/where/you/installed/libarpack.a"
before running Octave and try again

> B

c.

_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[Carlo de Falco-2]

On 16/ott/08, at 21:12, Ilan Boaz wrote:

>
> On Oct 16, 2008, at 10:06 AM, Carlo de Falco wrote:
>
>>
>> On 16/ott/08, at 17:29, Ilan Boaz wrote:
>>
>>>>
>>>> summarizing: you don't just copy manually the files from a  
>>>> package to the install location, you need to install them using  
>>>> the package manager pkg.m
>>>
>>> That was the first thing I tried, but it breaks -- see below. It  
>>> seems to be looking for eigs.oct, which was why I unpacked the  
>>> package to look for it.
>>
>> it tries to produce eigs.oct from the source code which is in the  
>> package and fails because it does not know where libarpack.a is as  
>> indicated by the error message:
>>
>>> "/usr/libexec/gcc/powerpc-apple-darwin8/4.0.1/ld: can't locate  
>>> file for: -larpack
>>> collect2: ld returned 1 exit status"
>>
>> so no eigs.oct is created
>>
>>> The thread
>>>>> https://www-old.cae.wisc.edu/pipermail/help-octave/2007-November.txt
>>> also suggests that eigs.oct should be inside this package.
>>
>> eigs.oct is not in the package tarball it will be in the package  
>> installation directory after compilation but you won't get it there
>> until you manage to build it
>>
>>> Any ideas?
>>
>> Set the environment variable LDFLAGS to point to
>> LDFLAGS="-Lthe/directory/where/you/installed/libarpack.a"
>> before running Octave and try again
>
> Thanks Carlo. That helped, but now:
> octave-3.0.2:1> pkg install  arpack-1.0.6.tar.gz
> /usr/libexec/gcc/powerpc-apple-darwin8/4.0.1/ld: Undefined symbols:  
> _G77_second_0
> collect2: ld returned 1 exit status
> make: *** [eigs.oct] Error 1
>
> ARPACK/README hints that the problem might be in the compilation of  
> libarpack.a:
> 6. You will also need to change the file "second.f" in the UTIL  
> directory
>   to whatever is appropriate for timing on your system.  The  
> "second" routine
>   provided works on most workstations.  If you are running on a Cray,
>   you can just edit the makefile in UTIL and take out the reference to
>  "second.o" to use the system second routine.
>
> I am guessing that ARPACK/UTIL/second.f does not work on my system,  
> hence the
> Undefined symbols: _G77_second_0
> I did not change ARPACK/UTIL/second.f, because I do not know how.
> I am using a Mac G5 PPC.
> I have g77 3.4.3. installed (using Fink Xcode 2.5), but there is no  
> second.f on my system.
> Web-searching for second.f for Mac yielded nothing.
>
> Any clue what to do?

I am not sure whether this is connected to your problem, but I usually  
try to avoid mixing libraries produced by different compilers, so as I  
am using the binary package for Octave 3.0.2 which is compiled with  
F77=fort77 I also compiled arpack with fort77.

Are you also using the binary package from octave-forge? In that case  
fort77 is included as
/Applications/Octave.app/Contents/Resources/bin/fort77

> B
c.

_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave