for loop with ode45 function

Hi,

I am trying to vary one of the parameters of my system. Then, I would like to input this value into the ode45 function, and determine the maximum displacement of my system, such that it is recorded in a variable, max_displ. I would like to iterate for many different values of my parameter, k_l, so I thought to use a for loop.

dk=1;
k_l=[1:dk:2];                                                               %0 to 2 is just a test 

for ii=1:dk:numel(k_l)                                                      %will iterate for all the values of the k_l from 0 up to 2
    k_l=k_l(ii)
    [t,z] = ode45(@(t, z) f(t,z,k_l,disp_data,vel_data,t_data),T,IC);
    disp_data_int=interp1(t_data,disp_data,t,'spline');                     %extracting the displacement data from the earthquake 
    displ=(z(:,8));                                                         %extracting the displacement data for a particular stiffness value 
    disp_a=abs(disp_data_int-displ);                                        %finding absolute value of diffence between displacment of earthquake and building 
    max_displ(ii)=max(disp_a)                                               %finding the maximumm displacement for a particular stiffness value 
end 
 
figure;
hold on 
plot(k_l,max_displ);

MATLAB indicates this error 'Index exceeds the number of array elements (1)', for this line of code,

k_l=k_l(ii)

Any suggestions would be greatly appreciated.

Thanks,

Kostas

Answer 1

[Star Strider]

Put the reference in the ‘f’ argument list:

[t,z] = ode45(@(t, z) f(t,z,k_l(ii),disp_data,vel_data,t_data),T,IC);

↑ ← HERE

and the loop indexing should be defined simply as:

for ii=1:numel(k_l)

.

[Kostas Kalabokas]

Thank you very much, that makes a lot more sense now.

K

[Star Strider]

As always, my pleasure!

[Kostas Kalabokas]

If I wanted to extend this to vary another parameter, lets say c_l, would I use a nested for loop, definining

for ii=1:numel(k_l)                                                         %will iterate for all the values of the k_l from 0 up to 2
    for jj=1:numel(c_l)
        k_lii=k_l(ii)
        c_ljj=c_l(jj)
        [t,z] = ode45(@(t, z) f(t,z,k_l(ii),c_l(jj),disp_data,vel_data,t_data),T,IC);
        disp_data_int=interp1(t_data,disp_data,t,'spline');                     %extracting the displacement data from the earthquake 
        displ=(z(:,8));                                                          %extracting the displacement data for a particular stiffness value 
        disp_a=abs(disp_data_int-displ);                                          %finding absolute value of diffence between displacment of earthquake and building 
        max_displ(ii,jj)=max(disp_a)                                                %finding the maximumm displacement for a particular stiffness value 
    end 
end 

[Star Strider]

Yes.

On a slightly different topic, I am not certain what you are doing with the interp1 call. If you want to have the results presented at specific time values, create ‘tspan’ (that appears to be ‘T’ here) as a vector of specific time instants, rather than a two-element range from minimum to maximum simulation times. The numerical solvers will then do that interpolation internally and present the results only at the indicated times.

[Kostas Kalabokas]

The interpolation function has been used to ensure that the two vectors t, and t_data are of the same size for the disp_a to be computed. It isn't a matter of specific time instants, but instead, the t_data has an interval of 0.01 seconds (over a 30 second time span). The interp1 fucntion is used to make sure that the vectors match.

If you think I can improve this further, then I would be very grateful of any recommendations. :)

K

[Torsten]

As said, if you specify T = t_data in the call to ode45, you don't need interp1. The returned t vector from the solver will match your t_data.

[Steven Lord]

Look at what happens when you specify a vector with more than 2 elements as the tspan input to ode45.

tspan = 0:0.25:5;
[t, y] = ode45(@(t, y) y, tspan, [0; 1]);
wereTspanValuesUsed = isequal(tspan.', t) % true

wereTspanValuesUsed = logical
   1

Alternately you could call ode45 with one output and use deval to evaluate the solution at the desired times.

sol = ode45(@(t, y) y, [0 5], [0; 1]);
norm(y - deval(sol, tspan).')

ans = 1.5993e-14

That's a pretty small difference.

Answer 2

[Torsten]

Look at what your loop does with the array k_l.

k_l is set to [1 2] before the loop.

For ii=1, the command

k_l = k_l(1)

sets k_l to a scalar, namely k_l(1) = 1.

For ii=1,

k_l = k_l(2)

tries to access the 2nd element from k_l.

But k_l is a scalar, thus has no second element. That's why MATLAB errors.

Workaround:

Replace

k_l = k_l(ii)

by

k_lii = k_l(ii)

[Kostas Kalabokas]

Thank you, that makes much more sense now.

Kostas