How can i pass an index to ode function?

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

Answer 1

[Star Strider]

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

[Francesco]

Hi Star, thank you for the quik answer. I did sm arrangement and it seems to work. Just 2 questions. Why you create A as A = randi(50,1,5)-1; and hy you put k1=1:size(A,2)?

thank you again Francescp

[Star Strider]

I defined A as I did because I wanted to be sure the code worked and the values of A were being incorporated into the function. The randi function output does not include zero, so to include zero as a possibility, I subtracted 1 from its output. Large integer random values of A will tell me more easily than small random values of A that the code works. I assume you have your own A vector, so you will substitute it for mine. Since A is a vector, I could have substituted length(A) for size(A,2) without altering the code functionality. The idea in that loop is that it chooses a value from A for each iteration and then integrates the DE. When it has run through all the values of A that are available, it stops and plots the results. The tv vector and yv matrix are the results. It was an interesting problem.

[Francesco]

Thank you so much for your help. I used it and work greatly!

[Star Strider]

My pleasure!

Here on MATLAB Answers, the sincerest form of thanks is to Accept the Answer that solves your problem. This also helps others who are searching here to find Answers that work for problems similar to theirs.

[Francesco]

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]

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.

Answer 2

[Sara]

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);

Answer 3

[Jan]

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