error using ODE45 to solve a 2nd order ODE
Show older comments
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
Accepted Answer
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
Tanya Shaban
on 6 Apr 2020
Star Strider
on 6 Apr 2020
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
on 7 Apr 2020
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
on 8 Apr 2020
Star Strider
on 8 Apr 2020
As always, my pleasure!
More Answers (0)
Categories
Find more on Data Type Identification in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
