Re: Help need for Differential algebraic equations

[Przemek Klosowski]

genehacker <address@hidden> tried to solve a differential equation of 
this form:

function res = f(x,xdot,t);
res = zeros(3,1);
k12=0.1;
k23=0.7;

xdot(1) = -k12*x(1);
xdot(2) = k12*x(1) - k23*x(2);
xdot(3) = k23*x(2);

endfunction;

This is not going to work, because you are forcing your xdot rather than
asking to find a value that satisfies your governing equation. Look at it
this way: you aren't using the return value res at all---you set it to
zero and the solver then can't progress with the solution and gets stuck
at zero.

Let's do a simple example: take a simple harmonic oscillator equation
of motion: x" + k*x=0, which can be transformed to a first-order by
substituting X=[x,x'], and rewritten as

                   
                Xdot + [ 0 k ; -1 0 ] * X = 0

To solve this with DASSL e.g. for k=1, you need to specify a
reasonable time interval that encompasses few oscillations, say from 1
to 10, and take an initial condition of zero speed and acceleration,
and non-zero initial 

     X0=[1 0]; Xdot0=[0 0]; t=0:.01:10;

and define the function in a way that returns the residual error
given the X and Xdot values that the solver might try:

 function res=f(X,Xdot,t);
 k=1; 
 res = Xdot + [0 k; -1 0] * X;
 endfunction

The solver must be called like this:

 [X,Xdot]=dassl('f',X0,Xdot0,t);

and indeed you see the harmonic motion with period 2pi: 

 plot(t/2/pi,X)

In your case, you probably need to rewrite your function as:

function res = f(x,xdot,t);

k = [-0.1  0     0 ; 
      0.1 -0.7   0 ; 
      0    0.7   0];

res=xdot-k*x;

endfunction