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)myodeError 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:
clcclearclose 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 onplot(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
Best Answer