Event function with odesolver error

3 次查看（过去 30 天）
Arun
Arun 2011-6-14
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.
  2 个评论
Matt Tearle
Matt Tearle 2011-6-14
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
Matt Tearle 2011-6-14
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?

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Matt Tearle
Matt Tearle 2011-6-14
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.
  4 个评论
显示 2更早的评论
Matt Tearle
Matt Tearle 2011-6-14
Syntactically that would be correct. (Can't tell you if it's mathematically/algorithmically correct!)

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by