Error when I run ODE45 function

Question

Satish Jawalageri 2020-7-15

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/565607-error-when-i-run-ode45-function

评论： Gifari Zulkarnaen 2020-7-27

Could I know where I might be wrong in ode.m as I am getting following message as error when I run the ODE file. F1 and F2 are function of time.

Error using odearguments (line 93)

OUTPUT must return a column vector.

Error in ode45 (line 115)

odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in Ode (line 11)

[t, x] = ode45(@Output,tspan,IC);

State space model: (Output.m)

function [sdot] = Output(t,x)
Ze = zeros(3);
I = eye(3);
M = [8070000,0,-629468070;0,8070000,112980;-629468070,112980,6.800000000000000e+10];
Ad = [8.053095218400001e+06,0,-4.831857131040000e+08;0,2.167940435676214e+05,0;-4.831857131040000e+08,0,3.865485704832000e+10];
C = [0,0,0;0,3.241885080000000e+05,0;0,0,1.301151158327999e+09];
Cm = [4.12e+04,0,-2.82e+06;0,1.19e+04,0;-2.82e6,0,3.11e+08];
MA = M+Ad;
Q = C+Cm;
Fwt = 2.793977550000000e+04;    %Wind force on turbine tower
Fg = -79086000;                 %Gravitational force
Fbuoy = 7.844814740000000e+07;  %Buoyancy force
Fp = 2.712318560000001e+06;     %Heave force
Fmx = -334731.8545;             
Fmz = -3517000;                 
t = 0:0.1:10
for i = 1:length(t)
    
    Fh = hydro(t(i))
    F1(1,i) = Fmx+27939.6+6.5*10^5+Fh(1,i);
    F2(1,i) = Fmz+Fg+Fbuoy+Fp;
    F3(1,i) = -112510430.2+3.44*10^6+266.5*10^5+(Fh(1,i)*18);
 A = [Ze, I; -inv(MA)*Q, Ze];
 B = [Ze; inv(MA)];
 F = [F1;F2;F3];
  
 Svec = [x(1);x(2);x(3);x(4);x(5);x(6)];
% 
 sdot = A*Svec + B*F;
end
end 

ODE.m:

%Initial Condition 
x(10) = 0; x(20) = 0; x(30) = 0;
x(40) = 0; x(50) = 0; x(60) = 0;
IC = [x(10);x(20);x(30);x(40);x(50);x(60)];
%time span
t0 = 0; tf = 10;
tspan = t0:0.1:tf;
[t, x] = ode45(@Output,tspan,IC);
IC = [IC;x(end)];
%x(1) = x(:,1);
%x(2) = x(:,2);
%x(3) = x(:,3);
%x(4) = x(:,4);
%x(5) = x(:,5);
%x(6) = x(:,6);
%Plot Results 
figure(1), clf 
plot(t,x(:,1)), xlabel('time(s)'), ylabel('surge(m)')
title ('Surge vs Time')

8 个评论
显示 6更早的评论隐藏 6更早的评论

Satish Jawalageri 2020-7-15

编辑：Satish Jawalageri 2020-7-15

在 MATLAB Online 中打开

Thanks for your reply. Here is the code:

hydro.m:

function [Fh] = hydro(t)
Cd = 0.6;
Ca = 0.9699;
R = 4.7;
t=0:0.1:10;
for i=1:length(t)
    vel = compute_wavevelocity(t(i));
    acc = compute_waveacceleration(t(i));
    fun = sym((0.5*997*Cd*2*R*abs(vel)*vel)+(Ca*997*pi*(R^2)*acc)+(pi*(R^2)*997*acc));
    Fh(1,i) = double(int(fun,0,-120));
    
end
plot(t,Fh), xlabel('time(s)'), ylabel('Fh')
title ('Fh vs Time')
end

compute_wavevelocity.m:

function velocity = compute_wavevelocity(t)
       T = 10;
       WN = (9.8*(T^2))/(2*pi());
       k = (2*pi())/WN;
       omega = (2*pi())/T;
       y = [320;310;300;290;280;270;260;250;240;230;220;210;200];
       d = 320;
       H = 4;
       vel = (omega*H/2)* ((cosh(k*y))/(sinh(k*d)))* cos((k*0)-(omega*t));
       
       velocity = mean(vel,1);
       
end

compute_waveacceleration.m:

function acceleration = compute_waveacceleration(t)
        T = 10;
       WN = (9.8*(T^2))/(2*pi());
       k = (2*pi())/WN;
       omega = (2*pi())/T;
       y = [320;310;300;290;280;270;260;250;240;230;220;210;200];
       d = 320;
       H = 4;
        acc = ((omega^2)*H/2)* ((cosh(k*y))/(sinh(k*d)))* sin((k*0)-(omega*t));
        
        acceleration = mean(acc,1);
               
    end

Walter Roberson 2020-7-17

Why are you looping over time there? why not just process for the one time that was passed into the function? The processing for any one t in the loop does not appear to depend on previous iterations.

Satish Jawalageri 2020-7-17

Till output.m file I have used looping for getting values corresponding the time. So, do I need to loop for ODE to process according to time?If so, how can I do?

Thanks.

请先登录，再进行评论。

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