MATLAB Answers

error using ODE45 to solve a 2nd order ODE

1 view (last 30 days)
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

  0 Comments

Sign in to comment.

Accepted Answer

Star Strider
Star Strider on 5 Apr 2020
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.

  6 Comments

Show 3 older comments
Star Strider
Star Strider on 7 Apr 2020
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
Tanya Shaban on 8 Apr 2020
That worked! Thank you so much for all your help!

Sign in to comment.

More Answers (0)

Sign in to answer this question.

Tags