How can i pass an index to ode function?

5 次查看(过去 30 天)
Hi all,
I'm a new user and I'm having trouble with ode45 function for a second order differential equation.
I should find a soultion for the problem d2y/dt2=(1-A*y)dy/dt-y in which A is a vector that depends on the time. I use [t,f]=ode45('fun',t,y0); How can I tell to the software that at timestep 1 I should use A(1) or a timestep 2 I should use A(2)?
thank you in advance Francesco
thank u advance

采纳的回答

Star Strider
Star Strider 2014-4-24
I defined ‘A’ as a large-amplitude random vector to be sure it worked:
A = randi(50,1,5)-1; % Create ‘A’
fun = @(t,y,A) [y(2); (1 - y(1) - A.*y(2))];
ti = 0:0.01:0.09; % Incremental time vector
te = 0; % Initial value for ‘t(end)’
ye = eps*[1 1]; % Initial ‘initial conditions’ vector
tv = []; % Initialise ‘tv’ to accumulate ‘t’
yv = []; % Initialise ‘yv’ to accumulate ‘y’
for k1 = 1:size(A,2)
t = te + ti; % Increment ‘t’ for this iteration
a = A(k1); % Select new value for ‘A’
[t, y] = ode45(@(t,y)fun(t,y,a), t, ye);
te = t(end);
ye = y(end,:); % Becomes new initial conditions
tv = [tv; t]; % Add current run to previous ‘tv’ vector
yv = [yv; y]; % Add current run to previous ‘yv’ matrix
end
figure(1)
plot(tv, yv)
legend('Y_1', 'Y_2', 'Location', 'NW')
xlabel('T')
grid
  6 个评论
Francesco
Francesco 2014-5-7
Hi, sorry if I disturb you again but I complicated my script since I should solve a system of equation. Following your suggestion I wrote this script (same stuff are not really useful but I put it anyway)
fun = @(tx,x,A,Bx,Dx,Cx,Ex,F) [x(2); -A.*((Bx.*c_d.*U_eff+Cx)*x(1)-(Dx.*c_l+Ex-F*cos(teta)))];
fun1 = @(ty,y,A,By, Dy, Cy,F) [y(2); A.*(By.*c_d.*U_eff.*(U_l-y(2))-Cy*y(2)+Dy-F*sin(teta))];
ti = 0:1:8; % Incremental time vector te = 0; % Initial value for ‘t(end)’ xe = [0 0]; % Initial ‘initial conditions’ vector ye = [0 0]; % Initial ‘initial conditions’ vector tvx = []; % Initialise ‘tv’ to accumulate ‘tx’ tvy = []; % Initialise ‘tv’ to accumulate ‘ty’ yv = []; % Initialise ‘yv’ to accumulate ‘y’ xv = []; % Initialise ‘xv’ to accumulate ‘x’
teta=0;
for k1 = 1:length(A)
tx = ti+te; % Increment ‘t’ for this iteration
ty = ti+te;
a = A(k1); % Select new value for ‘A’
bx = Bx(k1);
by = By(k1);
dx = Dx(k1);
dy = Dy(k1);
cx = Cx(k1);
cy = Cy(k1);
ex = Ex(k1);
f = F(k1);
U_eff = sqrt((U_l-ye(2))^2+xe(2)^2)
[tx, x] = ode45(@(tx,x)fun(tx,x,a,bx,dx,cx,ex,f), tx, xe);
[ty, y] = ode45(@(ty,y)fun1(ty,y,a,by,dy,cy,f), ty, ye);
teta = y(1)/x(1)
te = t(end);
xe = x(end,:) % Becomes new initial conditions
ye = y(end,:) % Becomes new initial conditions
tvx = [tvx; tx] % Add current run to previous ‘tv’ vector
tvy = [tvy; ty] % Add current run to previous ‘tv’ vector
yv = [yv; y] % Add current run to previous ‘yv’ matrix
xv = [xv; x] % Add current run to previous ‘yv’ matrix
end
All the variable are defined previously.
When I run the code the error message is "??? Error using ==> odearguments at 117 Solving @(TX,X)FUN(TX,X,A,BX,DX,CX,EX,F) requires an initial condition vector of length 2." I checked the vector xe and its length is 2. Do you know where is the problem?
thank you in advance Francesco
Star Strider
Star Strider 2014-5-7
I can’t run your code because there are too many undefined variables (from A to Ex). Since the first call to ode45 seems to be throwing the error, I suggest you insert:
IC_xe = xe % Check initial condition vector ‘xe’
just before your initial call to ode45 (the one that calls fun). That will write the value of xe to the Command Window, and tell you what xe is.
It might be necessary to insert a similar line before the second call (the one that calls fun1) in case that throws a similar error. That way you’ll know what is being passed to ode45 in the following line.

请先登录,再进行评论。

更多回答(2 个)

Sara
Sara 2014-4-24
You will need check the time into the fun and pass both A and t to it. You want something like:
[t,f]=ode45(@(t,x)fun(t,x,t,A),t,y0);
function dy = fun(t,x,timeserie,A)
k = find(t >= timeserie,1);
k = max(1,k-1);
a = A(k);

Jan
Jan 2014-4-25
It is important, that the function to be integrated is smooth. Otherwise the result of the integration is suspicious. See http://www.mathworks.com/matlabcentral/answers/59582#answer_72047

类别

Help CenterFile Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by