Using ode45 within a for loop

Hi,

So I am trying to use ode45 function when a parameter is changing. I have tried to use a for loop, but MATLAB is displaying some errors, and I do not think I know how to set up the for loop properly.

The parameter I want to vary is c_l

syms c_l
C_n=[ 0  0  0  0  0  0  0  0  0   0      0;
    0  0  0  0  0  0  0  0  0   0      0;
    0  0  0  0  0  0  0  0  0   0      0;
    0  0  0  0  0  0  0  0  0   0      0;
    0  0  0  0  0  0  0  0  0   0      0;
    0  0  0  0  0  0  0  0  0   0      0;
    0  0  0  0  0  0  0  0  0   0      0;
    0  0  0  0  0  0  0  0  0   0      0;
    0  0  0  0  0  0  0  0  0   0      0;
    0  0  0  0  0  0  0  0  0  c_l   -c_l;
    0  0  0  0  0  0  0  0  0  -c_l   c_l];  
C=zeros(11);
C=C_b+C_n;

Such that I have defined the ode45 function in the following way,

dz=@(t,z,c_l) [                                                   z(12:22);
                inv(M)*( k*disp(t)+C(1,1)*vel(t)-[0;0;0;0;0;0;0;0;0;knl*(z(10)-z(11))^3;(-knl)*(z(10)-z(11))^3]-C*z(12:22)-K*z(1:11))];
T=[0 30];
IC=zeros(22,1);

I am aware that the variable c_l does not appear directly in the dz function, but it is in the C_n matrix, is this the problem?

When a specific value of c_l is chosen, the code works fine. However, I would like to plot the function for varying values of c_l, and have tried this method

for c_l = 0:1*10^4:5*10^4;
    [t,z] = ode45(dz,T,IC);
    figure;
    plot(t, z(:,1))
    xlabel('t [s]')                    
    ylabel('x [m]') 
end 

Any suggestions would be greatly appreciated!

Answer 1

[Steven Lord]

While you could solve the problem using a complicated anonymous function and symbolic variables, I'd recommend instead writing a function file and skipping the (very simple) symbolic calculations. This will let you split the complicated definition of dz into smaller pieces that are easier to explain with comments and hopefully easier to debug if necessary.

T=[0 30];
IC=zeros(22,1);
for c_l = 0:1*10^4:5*10^4;
    [t,z] = ode45(@(t, z) f(t, z, c_l),T,IC);
    figure;
    plot(t, z(:,1))
    xlabel('t [s]')                    
    ylabel('x [m]') 
end 

function dz = f(t, z, c_l)
C = zeros(11);
C(10:11, 10:11) = c_l*[1 -1; -1 1];
% Add code here to define C_b or specify it as a fourth input
C = C_b + C;
% Define M, k, vel, knl, and K here
%
% Consider not using both lower-case k and upper-case K in your code
%
% Don't define a variable named disp or overload that name with
% a function that does calculations. It already has a meaning in MATLAB.
% I removed it from the line defining Q, replacing disp(t) with just t.
vec = [zeros(9, 1);
       knl*(z(10)-z(11))^3;
       (-knl)*(z(10)-z(11))^3];
Q = k*t+C(1,1)*vel(t)-vec-C*z(12:22)-K*z(1:11);
dz = [z(12:22);
      M\Q];
end

[Kostas Kalabokas]

Hi,

Thanks for replying! There is a lot of code defininig the M, K, and C matrices, so I have put this all into a function. I have extracted data from a frequency spectrum to define the displacement (displ(t)) and velocity (vel(t)), which I have also had to include in the function.

function dz=f(t,z,c_l) ... 
    (%contains all the code)
end

I have also defined dz in the way you recommended. Although there are error messages popping up that it does not work, as there are not enough input arguments.

T=[0 30];
IC=zeros(22,1);
for c_l = 0:1*10^4:5*10^4;
    [t,z] = ode45(@(t, z) f(t, z, c_l),T,IC);
    figure;
    plot(t, z(:,1))
    xlabel('t [s]')                    
    ylabel('x [m]') 
end 

I have then used the for loop that you defined in your answer.

When specifying a function does the function need to be saved into another file, if so, how do I call it into another script?

In addition, the for loop still does not work.

Thanks,

Kostas