error using ODE45 to solve a 2nd order ODE

Hi, I am not very experienced with MATLAB and im trying to solve a very long ODE. I've been playing with it for a while now trying to get it to work, and I cant seem to get past this error. The code is particularily ugly and complicated so I apologise if it is hard to navigate- I'd appreciate any suggestions! The error that I am receiving looks like this:

Not enough input arguments.

Error in flappingflightmodel>myode (line 102)
dpsidt(1) = psi(2);

Error in flappingflightmodel>@(t,psi)myode

Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

and my entire code is:

clc
clear
close all
%% Inputs
% rhat 0.4 AR 2
R = 0.3;            % wing lenth in m
x_LE_hat = 0.5;   % non-dimensinal position of the wing leading edge
yr = 0.0;           % y dimension of wing root offset in m (along span)
xr = 0.1;           % x dimension of wing root offset in m (along chord)
AR = 2;           % single wing aspect ratio
r1_hat = 0.4;    % non-dimensinal radius of first moment of wing area
f=7; %hz% waveform frequency 
T=1/f;% time period of one flapping cycle
n=1000; %points to be taken
dt=T/n;
phimax=pi/2;%maximum flapping angle
%the piching angle (theta)is 0 assuming no out of plane motion
psimax= pi/4;% max rotation angle
%AIR PARAMETERScurrent workspace and function context to the workspace and function context of the calling function or script in debug mode. 
rhoa=1.225; %air density at sea level (kg/m^3)
v=18.03E-6;%KINEMATIC VISCOSITY
%% Beta function
c_bar=R/AR;                     % mean geometric chord
yr_hat=yr/c_bar; %non-dimensional y-root offset
xrhat=xr/R; %non-dimensional x-root offset
n = 1000;                         % no of points to draw
r2_hat = 0.929*(r1_hat^0.732);    % from Ellington statistical relation
p = r1_hat*(((r1_hat*(1-r1_hat))/((r2_hat^2)-(r1_hat^2)))-1);
q = (1-r1_hat)*(((r1_hat*(1-r1_hat))/((r2_hat^2)-(r1_hat^2)))-1);
F = @(x)x.^(p-1).*((1-x).^(q-1));
B = quad(F,0,1);
dr = R/n;
for i=1:(n+1)
    r(i)=(i-1)*dr;
    r_hat(i)=r(i)/R;
    c_hat1(i)=(((r_hat(i))^(p-1))*((1-(r_hat(i)))^(q-1)))/B;
    c(i)=c_hat1(i)*c_bar;
    x_LE(i)=yr+r(i);
    x_TE(i)=yr+r(i);
    y_LE(i)=xr+x_LE_hat*c(i);
    y_TE(i)=xr+(x_LE_hat-1)*c(i);
end
%% Plotting
% subplot(3, 3, 1)
plot(x_LE,y_LE,'-','LineWidth',2,'Color','red','MarkerEdgeColor','red','MarkerFaceColor','red','MarkerSize', 10)
hold on
plot(x_TE,y_TE,'-','LineWidth',2,'Color','red','MarkerEdgeColor','red','MarkerFaceColor','red','MarkerSize', 10)
title('wing planform shape')
yhat_LE=y_LE/c_bar; %leading edge profile
%inertia tensor calc
    for i=1:n+1
        y(i)=y_LE(i)-y_TE(i);
    end
     y2=y.^3; %
     Ixx=sum(y2.*dr);
     Iyx=sum(y.^2.*r.*0.5.*dr);
     Ixy=Iyx;
% Iyy and Izz components get NEGLECTED---unneeded for the purposes of this model
t=linspace(0,T,n);
%kinematics
phi=phimax*cos(2*pi*f*t);%instantaneous flapping angle 
phidot=-2*pi*sin(2*pi*f*t);
phidotdot=-4*pi^2*phimax.*cos(2*pi*f.*t);
ubar=2*2*phi*f*R;%wing mean translational velocity
Re=(ubar*c_bar)/v;%reynolds number for flapping flight
%HINGE PARAMETERS
Eh=0.76E6; %wing membrane polymer stiffness
th=75E-6; %polymer thickness 
wh=R;%hinge width
Lh=0.1*(max(y_TE-y_LE));%flexural hinge length (assume 10% of wing max chord)
omegah=-phidot; %hingeline angular velocity assuming no stroke plane deviation
kh=(Eh*th^3*wh)/12*Lh;%hinge polymer stiffness
%%AERODYNAMIC FORCES
Fhat=(r2_hat^2)+(2*xrhat*r1_hat)+(xrhat^2);%non dimensional aerodynamic force 
yh_hat=0.5*c_hat1-yhat_LE-yr_hat; %non-dimensional rotational axis offset
Ixy_am=sum((r_hat+xrhat).*c_hat1.^2.*yh_hat);
Ixx_am=sum(c_hat1.^2.*(yh_hat.^2+(1/32).*c_hat1.^2));
%% SOLVING THE EQUATION OF MOTION
%solving through runge kutta methods
tspan = [0 T];
y0 = [2*pi 0.01];
[t,psi] = ode45(@(t,psi)myode, tspan ,y0);
plot(t,psi(:,1),'-o',t,psi(:,2),'-.')

function dpsidt = myode(t,psi,phidot,phi,phidotdot,Ixx,Ixy,c_bar,yr_hat,Fhat,rhoa,c_hat1,yhat_LE,Ixx_am,Ixy_am,R,f)
%the value of psi(1) is the instantaneous value of the rotational angle and psi(2) is the first derivitive 
Clmax=1.8;
CD0=0.4;
CDmax=3.4;
C_rd=CDmax;%rotational damping moment coefficient (should range from 3-6 for best agreement between measured and predicted results)
dpsidt = zeros(2,1);
dpsidt(1) = psi(2);
dpsidt(2) =(-(pi/4*rhoa*c_bar^3*R^2*Ixy_am.*(-phidotdot).*cos(psi(1))))+-0.5*rhoa.*psi(2).*abs(psi(2))*C_rd*c_bar^4*R(sum(0.25*((abs(yr_hat+yhat_LE).*((yr_hat+yhat_LE).^3))-((abs(yr_hat+yhat_LE-c_hat1)*((yr_hat+yhat_LE-c_hat1).^3))))))+sign((pi/2)-psi(1))*0.5*rhoa.*-phidot.^2*((Clmax*sin(2*(pi/2)-psi(1)))*cos((pi/2)-psi(1))+((CDmax+CD0)/2)-(((CDmax-CD0)/2)*cos(2*(pi/2)-psi(1)))*sin((pi/2)-psi(1)))*c_bar^2*R^3*Fhat*(sum((yr_hat+y_LE-c_hat1*((1/pi)*0.82*abs(psi(1))+0.05))*(r_hat+xrhat)^2.*c_hat1))+Ixy.*phidotdot.*cos(psi(1))+0.5*Ixx.*phidot^2*sin(2.*psi(1))/(Ixx+(pi/4*rhoa*c_bar^4*R*Ixx_am));
end

Answer 1

[Star Strider]

Note that psi is the polygamma function. The code ‘overshadows’ that function with the variable name.

The actual problem is that the ode45 call needs to be:

[t,psi] = ode45(@(t,psi)myode(t,psi,phidot,phi,phidotdot,Ixx,Ixy,c_bar,yr_hat,Fhat,rhoa,c_hat1,yhat_LE,Ixx_am,Ixy_am,R,f), tspan ,y0);

so that all the arguments are present in the argument list in the ‘myode’ call.

The (partially corrected) code is now:

[t,psi] = ode45(@(t,psi)myode(t,psi,phidot,phi,phidotdot,Ixx,Ixy,c_bar,yr_hat,Fhat,rhoa,c_hat1,yhat_LE,Ixx_am,Ixy_am,R,f), tspan ,y0);
plot(t,psi(:,1),'-o',t,psi(:,2),'-.')
function dpsidt = myode(t,psi,phidot,phi,phidotdot,Ixx,Ixy,c_bar,yr_hat,Fhat,rhoa,c_hat1,yhat_LE,Ixx_am,Ixy_am,R,f)
%the value of psi(1) is the instantaneous value of the rotational angle and psi(2) is the first derivitive 
Clmax=1.8;
CD0=0.4;
CDmax=3.4;
C_rd=CDmax;%rotational damping moment coefficient (should range from 3-6 for best agreement between measured and predicted results)
dpsidt = zeros(2,1);
dpsidt(1) = psi(2);
dpsidt(2) =(-(pi./4.*rhoa.*c_bar.^3.*R.^2.*Ixy_am.*(-phidotdot).*cos(psi(1))))+-0.5.*rhoa.*psi(2).*abs(psi(2)).*C_rd.*c_bar.^4.*R.*(sum(0.25.*((abs(yr_hat+yhat_LE).*((yr_hat+yhat_LE).^3))-((abs(yr_hat+yhat_LE-c_hat1).*((yr_hat+yhat_LE-c_hat1).^3))))))+sign((pi./2)-psi(1)).*0.5.*rhoa.*-phidot.^2.*((Clmax.*sin(2.*(pi./2)-psi(1))).*cos((pi./2)-psi(1))+((CDmax+CD0)./2)-(((CDmax-CD0)./2).*cos(2.*(pi./2)-psi(1))).*sin((pi./2)-psi(1))).*c_bar.^2.*R.^3.*Fhat.*(sum((yr_hat+y_LE-c_hat1.*((1./pi).*0.82.*abs(psi(1))+0.05)).*(r_hat+xrhat).^2.*c_hat1))+Ixy.*phidotdot.*cos(psi(1))+0.5.*Ixx.*phidot.^2.*sin(2.*psi(1))./(Ixx+(pi./4.*rhoa.*c_bar.^4.*R.*Ixx_am));
end

With those correctrions, the ‘dpsidt(2)’ assignment now throws:

Unable to perform assignment because the left and right sides have a different number of
elements.

You need to sort that, since it is not obvious from your code what that assignment is supposed to return. That line evaluates to a (1x1000) vector, which is as close as I can get to helping you sort it.

[Tanya Shaban]

Thank you!

I made that adjustment to the ode45 call and made the size of the arrays the same. When I run it now i get the error

Array indices must be positive integers or logical values
Error in flappingflightmodel>myode (line 104)
dpsidt(2)
=(-(pi/4*rhoa*(c_bar^3)*R^2*Ixy_am.*-phidotdot.*cos(psio(1)))+-0.5*rhoa*psio(2)*abs(psio(2))*C_rd*(c_bar^4)*R(sum(0.25*((abs(yr_hat+yhat_LE).*((yr_hat+yhat_LE).^3))-((abs(yr_hat+yhat_LE-c_hat1).*((yr_hat+yhat_LE-c_hat1).^3))))))+sign((pi/2)-psio(1))*0.5*rhoa*-1*(phidot.^2)*((Clmax.*sin(pi-psio(1))).*cos((pi/2)-psio(1))+((CDmax+CD0)/2)-(((CDmax-CD0)/2).*cos(2*(pi/2)-psio(1))).*sin((pi/2)-psio(1)))*(c_bar^2)*(R^3)*Fhat*(sum((yr_hat+y_LE-c_hat1*((1/pi)*0.82.*abs(psio(1))+0.05))*((r_hat+xrhat)^2)*c_hat1))+Ixy*phidotdot.*cos(psio(1))+0.5*Ixx*(phidot^2).*sin(2*psio(1)))/(Ixx+(pi/4*rhoa*(c_bar^4)*R*Ixx_am));

I'm not really sure what's going wrong. how can I move forward from here?

[Star Strider]

My pleasure!

This is the problem:

R(sum(0.25*((abs(yr_hat+yhat_LE)

Note that ‘R’ is a constant, so MATLAB interprets that expression as indexing into ‘R’. Assuming that you intend to multiply the expression in the parentheses by ‘R’, it should be:

R.*(sum(0.25*((abs(yr_hat+yhat_LE)

See if that (in addition to your other changes) now allows your function to integrate.

[Tanya Shaban]

Oh thanks! I completely missed that when I was checking through myself. I feel like I'm a lot closer to solving the problem now. I am still getting this prompt:

Unable to perform assignment because the left and right sides have a different number of elements.
Error in flappingflightmodel>myode (line 103)
dpsidt(2)
=(-(pi/4.*rhoa.*(c_bar^3).*R^2.*Ixy_am.*-phidotdot.*cos(psio(1)))+-0.5.*rhoa.*psio(2).*abs(psio(2)).*C_rd.*(c_bar^4).*R.*(sum(0.25.*((abs(yr_hat+yhat_LE).*((yr_hat+yhat_LE).^3))-((abs(yr_hat+yhat_LE-c_hat1).*((yr_hat+yhat_LE-c_hat1).^3))))))+sign((pi/2)-psio(1)).*0.5.*rhoa.*-1.*(phidot.^2).*((Clmax.*sin(pi-psio(1))).*cos((pi/2)-psio(1))+((CDmax+CD0)/2)-(((CDmax-CD0)/2).*cos(2*(pi/2)-psio(1))).*sin((pi/2)-psio(1))).*(c_bar^2).*(R^3).*Fhat.*(sum((yr_hat+y_LE-c_hat1.*((1/pi).*0.82.*abs(psio(1))+0.05)).*((r_hat+xrhat).^2).*c_hat1))+Ixy.*phidotdot.*cos(psio(1))+0.5*Ixx.*(phidot.^2).*sin(2*psio(1)))/(Ixx+(pi/4*rhoa*(c_bar^4)*R*Ixx_am));

the script currently looks like this:

clc
clear
close all
%% Inputs
% rhat 0.4 AR 2
R = 0.3;            % wing lenth in m
x_LE_hat = 0.5;   % non-dimensinal position of the wing leading edge
yr = 0.0;           % y dimension of wing root offset in m (along span)
xr = 0.1;           % x dimension of wing root offset in m (along chord)
AR = 2;           % single wing aspect ratio
r1_hat = 0.4;    % non-dimensinal radius of first moment of wing area
f=7; %hz% waveform frequency 
T=1/f;% time period of one flapping cycle
n=1000; %points to be taken
dt=T/n;
phimax=pi/2;%maximum flapping angle
%the piching angle (theta)is 0 assuming no out of plane motion
psimax= pi/4;% max rotation angle
%AIR PARAMETERScurrent workspace and function context to the workspace and function context of the calling function or script in debug mode. 
rhoa=1.225; %air density at sea level (kg/m^3)
v=18.03E-6;%KINEMATIC VISCOSITY
%% Beta function
c_bar=R/AR;                     % mean geometric chord
yr_hat=yr/c_bar; %non-dimensional y-root offset
xrhat=xr/R; %non-dimensional x-root offset
n = 1000;                         % no of points to draw
r2_hat = 0.929*(r1_hat^0.732);    % from Ellington statistical relation
p = r1_hat*(((r1_hat*(1-r1_hat))/((r2_hat^2)-(r1_hat^2)))-1);
q = (1-r1_hat)*(((r1_hat*(1-r1_hat))/((r2_hat^2)-(r1_hat^2)))-1);
F = @(x)x.^(p-1).*((1-x).^(q-1));
B = quad(F,0,1);
dr = R/n;
for i=1:n
    r(i)=(i-1)*dr;
    r_hat(i)=r(i)/R;
    c_hat1(i)=(((r_hat(i))^(p-1))*((1-(r_hat(i)))^(q-1)))/B;
    c(i)=c_hat1(i)*c_bar;
    x_LE(i)=yr+r(i);
    x_TE(i)=yr+r(i);
    y_LE(i)=xr+x_LE_hat*c(i);
    y_TE(i)=xr+(x_LE_hat-1)*c(i);
end
%% Plotting
subplot(3, 3, 1)
plot(x_LE,y_LE,'-','LineWidth',2,'Color','red','MarkerEdgeColor','red','MarkerFaceColor','red','MarkerSize', 10)
hold on
plot(x_TE,y_TE,'-','LineWidth',2,'Color','red','MarkerEdgeColor','red','MarkerFaceColor','red','MarkerSize', 10)
title('wing planform shape')
yhat_LE=y_LE/c_bar; %leading edge profile
%inertia tensor
    for i=1:n
        y(i)=y_LE(i)-y_TE(i);
    end
     y2=y.^3; %
     Ixx=sum(y2.*dr);
     Iyx=sum(y.^2.*r.*0.5.*dr);
     Ixy=Iyx;
% Iyy and Izz components get NEGLECTED---unneeded for the purposes of this model
t=linspace(0,T,n);
%kinematics
phi=phimax*cos(2*pi*f*t);%instantaneous flapping angle 
phidot=-2*pi*sin(2*pi*f*t);
phidotdot=-4*pi^2*phimax.*cos(2*pi*f.*t);
ubar=2*2*phi*f*R;%wing mean translational velocity
Re=(ubar*c_bar)/v;%reynolds number for flapping flight
%HINGE PARAMETERS
Eh=0.76E6; %wing membrane polymer stiffness
th=75E-6; %polymer thickness 
wh=R;%hinge width
Lh=0.1*(max(y_TE-y_LE));%flexural hinge length (assume 10% of wing max chord)
omegah=-phidot; %hingeline angular velocity assuming no stroke plane deviation
kh=(Eh*th^3*wh)/12*Lh;%hinge polymer stiffness
Fhat=(r2_hat^2)+(2*xrhat*r1_hat)+(xrhat^2);%non dimensional aerodynamic force 
yh_hat=0.5*c_hat1-yhat_LE-yr_hat; %non-dimensional rotational axis offset
Ixy_am=sum((r_hat+xrhat).*c_hat1.^2.*yh_hat);
Ixx_am=sum(c_hat1.^2.*(yh_hat.^2+(1/32).*c_hat1.^2));
%% SOLVING THE EQUATION OF MOTION
%solving through runge kutta methods
tspan = [0 T];
y0 = [2*pi 0.01];
[psio,t] = ode45(@(t,psio)myode(t,psio,phidot,phi,phidotdot,Ixx,Ixy,c_bar,yr_hat,Fhat,rhoa,c_hat1,yhat_LE,Ixx_am,Ixy_am,R,f,y_LE,r_hat,xrhat), tspan ,y0);
plot(t,psio(:,1),'-o',t,psio(:,2),'-.')
% omega=sqrt(omegax^2+omegay^2+omegaz^2);
% 
% L=0.5*rhoa*omegah^2*CL*c_bar*R^3*Fhat; %lift generation with time
% plot(L,t)
% D=0.5*rhoa*omegah^2*CD*c_bar*R^3*Fhat; %drag variation
% N=0.5*rhoa*omegah^2*CN*c_bar*R^3*Fhat;%variation of normal force 
% ROTATIONAL AERODYNAMIC FORCE CONTRIBUTION TO LIFT
% W0dot=phidotdot.*cos(psi).*r;
% non-zero added mass coefficients
% lamdaz=pi*rhoa*(0.5*(y_LE-y_TE))^2;
% lamdazw=-pi*rhoa*(y_LE-y_TE)^2*yh;

function dpsidt = myode(~,psio,phidot,~,phidotdot,Ixx,Ixy,c_bar,yr_hat,Fhat,rhoa,c_hat1,yhat_LE,Ixx_am,Ixy_am,R,~,y_LE,r_hat,xrhat)
%the value of psi(1) is the instantaneous value of the rotational angle and psi(2) is the first derivitive 
Clmax=1.8;
CD0=0.4;
CDmax=3.4;
C_rd=CDmax;%rotational damping moment coefficient (should range from 3-6 for best agreement between measured and predicted results)
dpsidt = zeros(2,1);
dpsidt(1) = psio(2);
dpsidt(2) =(-(pi/4*rhoa*(c_bar^3)*R^2*Ixy_am.*-phidotdot.*cos(psio(1)))+-0.5*rhoa.*psio(2).*abs(psio(2))*C_rd*(c_bar^4)*R*(sum(0.25.*((abs(yr_hat+yhat_LE).*((yr_hat+yhat_LE).^3))-((abs(yr_hat+yhat_LE-c_hat1).*((yr_hat+yhat_LE-c_hat1).^3))))))+sign((pi/2)-psio(1))*0.5*rhoa*-1.*(phidot.^2)*((Clmax.*sin(pi-psio(1))).*cos((pi/2)-psio(1))+((CDmax+CD0)/2)-(((CDmax-CD0)/2).*cos(pi-psio(1))).*sin((pi/2)-psio(1)))*(c_bar^2)*(R^3)*Fhat.*(sum((yr_hat+y_LE-c_hat1*((1/pi)*0.82.*abs(psio(1))+0.05)).*((r_hat+xrhat).^2).*c_hat1))+Ixy.*phidotdot.*cos(psio(1))+0.5*Ixx.*(phidot.^2).*sin(2.*psio(1)))/(Ixx+(pi/4*rhoa*(c_bar^4)*R*Ixx_am));
end

the arrays are all 1x1000 so I'm not sure why this is the case.

[Star Strider]

The value assigned to ‘dpsidt(2)’ must be a scalar. One of the values used to calculate it is a (1x1000) element vector, and that is causing the problem. You need to find out which one it is, and choose the appropriate element of it to use in ‘dpsidt(2)’.

[Tanya Shaban]

That worked! Thank you so much for all your help!

[Star Strider]

As always, my pleasure!