Nested loops, Arrays and ODE's - Confusion

I have been having trouble with this code for quite some time now. The gist of it is, I have a sort of hysteresis loop I have to plot. This is time dependent and requires the solution of two different equation at different points of time.

I got the result for the first order part of it when I tried it seperately. When I combine the second order equation I am getting errors.

My idea is that say for example, first 5 seconds are to be solved by the first order eqn. The solution is stored in y1 variable. Then the next ten second is solved using second order eqn and the values are stored in variable y2. The next 6 seconds is again first order equation. The solution is to be appended to the variable y1. Then next set to y2 and so on till the time duration is over.

I tried nested loops and am confused as to how to proceed next.

My code is given below. Along with the 2 functions I used.

Right now the error is ---> ------------------------------------------ Error using odearguments (line 93) F returns a vector of length 31, but the length of initial conditions vector is 1. The vector returned by F and the initial conditions vector must have the same number of elements. --------------------------------- I would be grateful if anyone can help me out with this. Let me know if there is any further clarification.

-----------------------------

%Assuming the Airfoil undergoes a Pitching motion of

% theta = 10 + 10 sin(0.1*T)

% M = 0.3

clc;

clear all;

global A F G H I J K L W

vel=99; %Velocity in m/s

speedOfSound = 340;

M=vel/speedOfSound;

b=0.2; % Blade Semichord

ti =0:1:30; % Time in Seconds

t=ti.'; % time is made to a column vector

thetaD = 15*(1-M^2);

lambda = 0.17-0.13*M;

alpha = 0.53+0.25*(sqrt(1-M^2)-1);

s = (pi + 5*pi*( (( 1-M^2)^0.285)-1))*(pi/180);

k = (pi*0.5 + 1.96*pi*( sqrt(1-M^2)-1))*(pi/180);

sigma = ((2*pi)/sqrt(1-M^2))*(pi/180);

W = (0.1*vel)/b; % omega

theta0 = 10;

theta1 = 10;

theta = theta0+theta1*sin(W*t); %Pitch angle in degrees

deltaCz = zeros(length(t),1);

y1 = zeros(length(t),1);

y2 = zeros(length(t),1);

al= zeros(length(t),1);

rl= zeros(length(t),1);

El= zeros(length(t),1);

d= zeros(length(t),1);

Po = (0.100*(1-M^8))/sqrt(1-M^2);

P1 = 0.1*M^4;

bl = 0.7*(1-M);

hl = -0.5+(1.5-M)*M^2;

d1=0;

diffCz = Po;

Wo = theta*vel;

W1 = b*10*W*cos(W*t);

for n = 1:length(ti)

    if ( theta(n)<thetaD)

        deltaCz(n) = 0;

    else

        deltaCz(n) = (Po-P1)*(theta(n)-thetaD)- bl*((exp(hl*(theta(n)-thetaD)))-1);

    end

    al(n) = 0.30+0.20*(deltaCz(n)).^2;
    rl(n) = (0.20+0.20*(deltaCz(n)).^2).^2;
    El(n) = -0.05*(deltaCz(n)).^2;
    d(n) = d1* deltaCz(n);

end

% A to H is for the first oder diffrential equation

A = (lambda*vel)/b;

B = A*diffCz*vel;

C = A*sigma*b*theta1*W;

D = (alpha*diffCz +d)* vel*theta1*W;

E = alpha*sigma*b*theta1;

F = B*theta0 +B*theta1;

G = B*theta1-E;

H = C+D;

%I to is for the second order differential equation

I = al*(vel/b);

J = rl*(vel/b).^2;

K = J*vel*diffCz;

L = El*(vel/b)*theta1*W;

store1 =1;

store2 =1;

m = 1;

while (m<length(t))

    for n = store1:length(t)

        if(theta(n)>thetaD)

            store2=n;

            m=n;

            break

        end

    end

    t = [store1 store2];

    init = 0;

    [t,y1] = ode23(@f,t,init);

    for n = store2:length(t)

    if(theta(n)<thetaD)

        store1=n;

        m=n;

        break
    end

    end
end

Cz = (s*b*Wo + k*b*W1 + vel* y1 + vel*y2)/(vel^2);

plot(theta, Cz);

title('Pitching Motion');

xlabel('Theta (in degrees) -->');

ylabel('Cz');

----------------------

function r = f(t,y1)

global A F G H W

r = -A*y1+F*sin(W*t) + G*cos(W*t)+H;

-------------------------------- function r = f1(t,y2)

global I J K L W

r = zeros(2,1);

r(1) = y2(2);

r(2) = -(I*y2(2) + J*y2(1) + K + L*cos(W*t));

[Jan]

Please format your code properly (see the "? Help" link). When you mention, that you get "errors" (it should be one error only?!), post the complete error message. We cannot guess, where which error occurs.

[Jan]

The code is still very hard to read. A clean copy of the complete error message would be helpful also.

Answer 1

[Jan]

The part of the error message you have shown us says:

 F returns a vector of length 31, but the length of initial conditions vector is 1.

This means, that the function to be integrated replies a vector of length 31, but you have provided a scalar initial value only. Please use the debugger to find out more: Type this in the command line

dbstop if error

and start your program again. When it stops you can check the current values of the locally used variables.

The different global variables with one letter names are a really bad idea. Better provide parameters by anonymous functions.