Event function with odesolver error

Hello,

I am trying to integrate a vector E subject to a first order differential equation. I want the integration to stop when a certain threshold is reached. Therefore I have constructed the following functions:

d=@(s,E)dPdEfull(s,E,t,U0,H0,mu,initial,final)

s is my integration variable. E is a vector in time, with a predetermined number of components. My event function is

e=@(s,E)Pifevent(s,E,t,U0,H0,mu)

with Pifevent defining the value, isterminal, and direction. Now when I evaluate both functions separately with some E vector, they work fine. But when I try to solve by using

[s E]=ode45(d,[0 10],E0,options);

with

options=odeset('Events',e);

and E0 the initial E vector, I get an error message

Attempted to access E(2); index out of bounds because numel(E)=1.

I don't know why this happens, especially since when I try to integrate without the event function and just do

    [s E]=ode45(d,[0 10],E0);

it works fine. It seems that the event function is not reading E as a vector. Maybe it is trying to evaluate the event function before the first integration supplies a vector E? I don't know how to fix this.

Any help would be appreciated.

Edit: I'm posting more code as per Matt's request.

d(0,E0)

returns a column vector of size 101x1, I won't post the output since it's too long.

[a,b,c]=e(0,E0)

returns a=-.2923, b=1, and c=0. Does that help?

Edit 2:E0 is 1x101. The code for dPdEfull isn't that bad, I'll go ahead and post it and Pifevent:

function dEds = dPdEfull(s, E, t, U0, H0, mu, initial, final)
dEds=zeros([1 numel(E)]);
for l=1:length(t)
  m=min(t):t(2)-t(1):t(l);
  a=Schrodinger(t,U0,H0,mu,E);
  c=Schrodinger(m,U0,H0,mu,E);
  b=c';
  x=a(final,initial);
  Y=a*b*mu*c;
  y=Y(final,initial);
  dEds(l)=2*imag(x*y);
end
dEds=dEds(:);
end

and Pifevent:

function [value, isterminal, direction] = Pifevent(s,E,t,U0,H0,mu)
u=Schrodinger(t,U0,H0,mu,E(end,:));
v=u(2,3);
v=abs(v);
v=v^2;
value=v-0.99;
isterminal=1;
direction=0;
end

The function Schrodinger just takes in those arguments and returns a unitary matrix (3x3 in this example). By the way t is also a vector.

[Matt Tearle]

Is it possible for you to post more of the code? I can't reproduce this problem with a simple example (other than doing something really obviously wrong). How about a simple diagnostic of

>> d(0,E0)

>> [a,b,c] = e(0,E0)

[Matt Tearle]

So E0 is 101-by-1 also? I'm guessing that means the code for dPdEfull is too hideous to post? Not sure it will help, but can you post your code for Pifevent?

Answer 1

[Matt Tearle]

Ah. I think I have an idea now. Check the error message trace. I'm guessing the error actually occurs within Schrodinger. When ode45 passes E into the ODE function dPdEfull or the event function Pifevent, it passes just the current timestep's value, so E is a column vector (even though the E returned by ode45 is in the form of a matrix). Inside dPdEfull you call Schrodinger:

a=Schrodinger(t,U0,H0,mu,E);

and here E is a vector and all is well. But inside Pifevent you call

u=Schrodinger(t,U0,H0,mu,E(end,:));

Because E is a (column) vector, E(end,:) is the last row, which is a scalar. If Schrodinger expects a vector as the last argument, that could cause an error like you're getting.

[Arun]

Okay, this kind of advice is exactly what I need. I think you are right, the error traces it to Schrodinger inside Pifevent. It does expect a vector for the last argument. Is there some way I can pass the entire row to it?

[Arun]

Oh, should I remove the (end,:) and just write E there? That may be right...

[Matt Tearle]

Syntactically that would be correct. (Can't tell you if it's mathematically/algorithmically correct!)

[Arun]

Thank you for your help.