ODE23s, Differential equation, problem solving

1 次查看（过去 30 天）
Iro Malefaki
Iro Malefaki 2020-5-27
Hello everybody. I have a serious problem solving the diff. equation, expressed below. Any help;
function [kap, alf0, Pstar, Astar]= nofilter()
imu=sqrt(-1);
rho=1000;RR=0.3;
Jb=0.1; %diastato
CC=0.1;%diastato
chr=0.1;
cht=0.1;
U=0.5;
span=0.34/2;
alam1=90;
alam2=90;
chm=0.5*(chr+cht);
span3d=2*span;
area=2*(chm*span);
AR=(span3d^2)/area;
AR=10;
B=14.75; %adiastato
BB=B*rho*U*0.34*chr*chr*RR; %diastato
pic=0.25;
II=Jb/pi*rho*area*chr*chr*chr; %adiastato
C=CC/1/2*pi*rho*area*chr*U^2; %adiastato
a=span3d*RR^2/chr*area;
b=RR^2/chr^2;
e=RR/chr;
d=RR*span3d/area;
alf0=linspace(1, 60, 2);
kap=linspace(0.025,0.2,2);
for i=1:length(alf0)
for j=1:length(kap)
om(j)=pi*kap(j)*2*U/chr; % reduced velocity k, Huxham, k=pi*khuxham
T(j)=2*pi/om(j);
alf0r(i)=alf0(i)*pi/180;
theo(j)=besselh(1,2,kap(j))./(besselh(1,2,kap(j))+imu*besselh(0,2,kap(j))); %theo function, lift deficiency factor
Re=U*chr*10^6; %reynolds
CD0=0.075/(log10(Re)+2)^2; %skin friction
fprintf('Calculating for %d\n', i, j);
[t,y]=ode23s(@(t,y) nofilter(t,y, i, j), [0 , 100],[0, 0]);
[Astar(i,j), Pstar(i,j)] = processsignal(t,y);
Pstar(i,j) = 0.5*BB.*(Pstar(i,j))./0.5*rho* 2*RR.*(Astar(i,j)).*cos(((Astar(i,j)))).*U^3;
end
end
function dYdt=nofilter(t,Y, i, j)
thi(i,j)=alf0r(i).*cos(om(j)*t);
k1=8*II+2*a*cos(thi(i,j)).*(cos(Y(1))-Y(1)*sin(Y(1)));
dYdt=[Y(2);(1/k1)*((Y(2)).*(-2*B+theo(j).*(AR/AR+2).*b.*cos(thi(i,j)).*(-cos(Y(1)+Y(1).*sin(Y(1)))))+...
(Y(2).^2)*2*a*cos(thi(i,j))*(2*sin(Y(1))*Y(1)+Y(1)*cos(Y(1)))-Y(1)*C+...
theo(j)*(AR/AR+2)*e*(thi(i,j).*cos(thi(i,j))+(3/4-pic)*2*thi(i,j).*cos(thi(i,j)))+...
d*(thi(i,j)+(1/2-pic)*2*thi(i,j).*cos(thi(i,j)))-e*(1/pi)*sin(thi(i,j))*(CD0+1.6*(thi-atan(2*e*Y(2).*(cos(Y(1))-Y(1).*sin(Y(1))))).^2))];
end
end
function [Aout,Pstar] = processsignal(t,y)
% process signal
Y = y(:,1);
Ydot = y(:,2);
k = 2:length(Y)-1;
Kmax = find( Y(k)>Y(k-1) & Y(k)>Y(k+1) );
Kmax = Kmax+1; %correct index
Kmin = find( Y(k)<Y(k-1) & Y(k)<Y(k+1) );
Kmin = Kmin+1; %correct index
ti_max = 0.5*(t(Kmax(1:end-1))+t(Kmax(2:end)));
ti_min = 0.5*(t(Kmin(1:end-1))+t(Kmin(2:end)));
N = 5; % number of final cycles for averagingaw
Aout = 0.5*(mean(Y(Kmax(end-N,end))) - mean(Y(Kmin(end-N,end))));
subplot(2,1,1), plot(t,Y,t(Kmax),Y(Kmax),'ob',t(Kmin),Y(Kmin),'ob')
subplot(2,1,2), plot(t(Kmax),Y(Kmax),'o',t(Kmin),abs(Y(Kmin)),'o')
Jint = find(t>=t(Kmax(end-N)) & t<t(Kmax(end)));
Pstar = (1/4)*trapz(t(Jint),Ydot(Jint).^2); % integral,1/4
end

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Mathematics 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by