Solving coupled second order differential equations using ODE

Hi,

I have the following second order and first order differntial equations.

syms t wn n d33 Ft meff R Cs x(t) v(t) Y
Dx=diff(x);
D2x=diff(x,2);
Dv=diff(v);
Eq1=D2x==((-wn^2)*x)+(n*(wn^2)*d33*v)+(Ft/meff);
Eq2=Dv==((-1/R*Cs)*v)-((n*d33*meff*(wn^2))/Cs)*diff(x,t);
[V,Subs]=odeToVectorField(Eq1,Eq2);
ftotal= @(t,Y,wn,n,d33,Ft,meff,R,Cs)[-(Cs.*Y(1))./R-(d33.*meff.*n.*wn.^2.*Y(3))./Cs;Y(3);Ft./meff-wn.^2.*Y(2)+d33.*n.*wn.^2.*Y(1)];
ic=zeros(4,1);
tspan=[0,1];
[T,Y]=ode45(@(t,y) ftotal(t,Y,wn,n,d33,Ft,meff,R,Cs), tspan, ic);
plot(T,Y)
grid

So far I've got this, however when I run the code I get the following error.

Index exceeds the number of array elements (1).

Error in sym/subsref (line 902)
            R_tilde =
            builtin('subsref',L_tilde,Idx);

Error in
Draft_Code>@(t,Y,wn,n,d33,Ft,meff,R,Cs)[-(Cs.*Y(1))./R-(d33.*meff.*n.*wn.^2.*Y(3))./Cs;Y(3);Ft./meff-wn.^2.*Y(2)+d33.*n.*wn.^2.*Y(1)]
(line 62)
ftotal=
@(t,Y,wn,n,d33,Ft,meff,R,Cs)[-(Cs.*Y(1))./R-(d33.*meff.*n.*wn.^2.*Y(3))./Cs;Y(3);Ft./meff-wn.^2.*Y(2)+d33.*n.*wn.^2.*Y(1)];

Error in
Draft_Code>@(t,y)ftotal(t,Y,wn,n,d33,Ft,meff,R,Cs)
(line 65)
[T,Y]=ode45(@(t,y)
ftotal(t,Y,wn,n,d33,Ft,meff,R,Cs), tspan, ic);

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);

Error in Draft_Code (line 65)
[T,Y]=ode45(@(t,y)
ftotal(t,Y,wn,n,d33,Ft,meff,R,Cs), tspan, ic);

F(t) is an input force which varies over time. The time domain is predefined as 0:0.0001:0.025 for the input force.

From the code im aiming to get values for the voltage, v against time or v against x (the displacement).

Any help would be much appreciated.

Thanks

Answer 1

[Star Strider]

None of the variables you are passing to ‘ftotal’ have been defined numerically (at least in the posted code).

They need to be in order to use a numerical integration function such as ode45.

[George Jones]

Thanks for your response,

I've ordered them and included the variables after 'ftotal'.

This is the complete code, I'm still managing to get the same error messages. I'm a beginner when it comes to ode, apologies if I'm still doing something blantantly wrong.

f=466.7;
t=0:0.0001:0.025;
a=160*10^3;
y=(160*10^3)+(a*sin(2*pi*f*t));
d=0.0127;
A=(pi*(0.0127^2))/4; %area
F=y*A; %pressure*area to obtain force
syms wn n d33 Ft meff R Cs x(t) v(t) Y
Dx=diff(x);
D2x=diff(x,2);
Dv=diff(v);
Eq1=D2x==((-wn^2)*x)+(n*(wn^2)*d33*v)+(Ft/meff);
Eq2=Dv==((-1/R*Cs)*v)-((n*d33*meff*(wn^2))/Cs)*diff(x,t);
[V,Subs]=odeToVectorField(Eq1,Eq2);
ftotal=  @(t,Y,wn,n,d33,Ft,meff,R,Cs)[-(Cs.*Y(1))./R-(d33.*meff.*n.*wn.^2.*Y(3))./Cs;Y(3);Ft./meff-wn.^2.*Y(2)+d33.*n.*wn.^2.*Y(1)]
wn=sqrt(k/meff);%natural frequency of stack
n=250; %stack layers 
d33=6.35*(10^-4);%dialectric constant
Ft=F-ff; %force input (force - friction)
Ft(Ft<=ff)=ff; %limiting force input to exclude static friction
meff=7.08*(10^-3); %effective mass
R=20; %Load Resistance
Cs=2.6*(10^-6); %stack capacitance
ic=zeros(4,1);
tspan=[0,1];
[T,Y]=ode45(@(t,y) ftotal(t,Y,wn,n,d33,Ft,meff,R,Cs), tspan, ic);
plot(T,Y)
grid
    

[Star Strider]

My pleasure!

All the numeric variables ‘ftotal’ uses as extra parameters must be defined before the ode45 call. There are still several that appear to be missing.

Also, MATLAB is case-sensitive, so y~=Y. Change that —

[T,Y]=ode45(@(t,y) ftotal(t,Y,wn,n,d33,Ft,meff,R,Cs), tspan, ic);

as well to —

[T,Y]=ode45(@(t,y) ftotal(t,y,wn,n,d33,Ft,meff,R,Cs), tspan, ic);

.