Solving a coupled system of differential equations with varying orders.

Hi

I have a system of three equations. Two of them are second order differentials and one if a first order. I am unsure of how to develop the matrix for solving with ode45. Usually you would develop an equation for the second derivative however if this was the case with the first order equation you would lose all data. The system is seen in the picture.

Can anyone help me out?

Answer 1

[Torsten]

y1'=y2
y2'=(-b*y2-k*y1-mu*(y1-y3)+f0)/m
y3'=(Lambda*y5-mu*(y3-y1))/Lambda
y4'=y5
y5'=(-B*y5-mu*(y3-y1)-K*y4+fs)/M

where

y1=x,y2=x',y3=y,y4=z,y5=z'

Best wishes

Torsten.

[York Stanham]

Hi Torsten.

Thanks for that solution. I had something similar however the expression for y5' has has the derivative of (z-y).

Would it be correct to write it as:

y5'=(-B*y5-lambda*(y5-y3')-K*y4+fs)/M

I have included a derivative (y3') on the right side of the expression. I am not sure this is correct.

Thanks, York

[Torsten]

It's correct, but ODE45 does not accept the right-hand side of equation 5 as written.

You have to take y3' to the left-Hand side and define a mass Matrix as

M=[1 0 0 0 0; 0 1 0 0 0; 0 0 1 0 0; 0 0 0 1 0; 0 0 -Lambda/M 0 1]

and the right-hand side as

dydt=[ y2;(-b*y2-k*y1-mu*(y1-y3)+f0)/m;(Lambda*y5-mu*(y3-y1))/Lambda;y5;(-B*y5-Lambda*y5-K*y4+fs)/M]

Best wishes

Torsten.

[York Stanham]

Thanks Torsten

I am still having a few issues. Running your code results in a matrix dimension error when using "/"

Any idea how I can fix this?

Thanks York

[Torsten]

Please show the MATLAB code you are using.

Best wishes

Torsten.

[York Stanham]

No problem. Here it is:

%define tstart, tend and the number of time steps
tstart=0;
tend=500;
n=1000;
tspan=linspace(tstart,tend,n);

%define initial conditions
xinit=[0;0;0;0;0];

%get x
[t,x]=ode45(@OWCsolver,tspan,xinit)
plot(t,x(:,1),t,x(:,3),t,x(:,4))

* _ | |Function| | _ *

function dxdt=OWCsolver(t,x)

%---------------------------- ODE45 Function -----------------------------%

%Constants
rhow=1025;                                   %Density of sea water (kg/m^3)
g=9.81;                                %Acceleration due to gravity (m/s^2)

%------------------------------- Structure -------------------------------%

%------------------------- Structure Dimensions --------------------------%

h=50;                                                 %Structure height (m)
wt=1;                                                   %Wall thickness (m)
Dia=60;                                             %Structure diameter (m)
R=h/2;                                          %Structure outer radius (m)
r=R-wt;                                         %Structure inner radius (m)
rhos=600;                            %Average density of structure (kg/m^3)

%------------------------------ Calculations -----------------------------%

areas=pi*(R^2-r^2);                           %Base area of Structure (m^2)
areao=pi*r^2;                                       %Base area of OWC (m^2)
draft=h*rhos/rhow;                                     %Structure draft (m)

m=areao*(draft)*rhow;                                        %OWC Mass (kg)
k=rhow*g*areao;                                        %OWC Stiffness (N/m)
b=0.02*sqrt(k*m);                             %OWC Radiation Damping (Ns/m)

mu=k/0.1;                                        %Air Compressibility (N/m)
lambda=1.225*areao^3/(2*0.64^2*areao^2);      %Power takeoff Damping (Ns/m)

M=h*areas*rhos;                                        %Structure Mass (kg)
K=0.01*mu*1000;                                  %Structure Stiffness (N/m)
B=0.02*sqrt(K*M);                                 %Structure Damping (Ns/m)

%--------------------------------- Wave ----------------------------------%

%-------------------------------- Inputs ---------------------------------%
L= 100;                                                    %Wave Length (m)
T= 10.8529;                                                %Wave period (s)
d= 100;                                                    %Water depth (m)
a= 1;                                                   %Wave amplitude (m)

%------------------------------ Calculations -----------------------------%

WN=2*pi/L;                                              %Wave number (m^-1)
omega=2*pi/T;                          %Wave radian frequency (radians/sec)
H=2*a;                                                     %Wave height (m)
y=a*sin(omega*t);                                   %Wave displacement (m)                
w=sqrt(k/m);                                %OWC Undamped Natural Frequency  
W=sqrt(K/M);                          %Structure undamped Natural Frequency

%-------------------------------- Matrix ---------------------------------%

%----------------------------- Define Matrix -----------------------------%
dxdt=zeros(5,1);   

%--------------------- Differential equation matrix ----------------------%
f0=areao/(areao+areas)*a*sin(omega*t);
fs=areas/(areao+areas)*a*sin(omega*t);

M=[1 0 0 0 0
   0 1 0 0 0
   0 0 1 0 0 
   0 0 0 1 0
   0 0 -lambda/M 0 1];

dxdt=[x(2);(-b*x(2)-k*x(1)-mu*(x(1)-x(3))+f0)/m;(lambda*x(5)-mu*(x(3)-x(1)))/lambda;x(5);(-B*x(5)-lambda*x(5)-K*x(4)+fs)/M];

end

[Torsten]

You have to specify M before you call ODE45 via odeset.

Take a look at the example hb1dae under

http://de.mathworks.com/help/matlab/math/ordinary-differential-equations.html#f1-670396

Best wishes

Torsten.