I need to to solve this coupled ODEs. I have written the code but not really getting any output. Any help would be appreciated!

function [time, x, xdot, y] = bluff

meff = 0.0114; %Effective mass (kg)
ceff = 0.0212; %Effective damping (N/(m/s))
keff = 8.86;   %Effective stiffnes (N/m)
thet = 0.000019758; %Electromechanical coefficient (N/V)
Cp = 0.0000000157; %Capacitance
rho = 1.225; %air density (kg/m3)
u = 10; %wind velocity(m/s)
stip = 0.0118; %exposed area of the bluff body(m2)
L = 0.15 ; %length of the beam
a1 = 2.3;  %empirical coefficient for the aerodynamic force calculation
R = 1050000; %load resistance (ohm)


%values
x0 = 1; xdot0 = 0; y0 = 1;
t0 = 0; tf = 60;

%Time span
tspan = [t0, tf];


%initial conditions
IC = [x0, xdot0,y0];

% [sdot] = g(t,s)
sdot = @(t, s) [s(2);
    (-ceff*s(2) - keff*s(1) - thet*s(3) + 0.5*rho*stip*(a1*((s(2)/u)+(1.5*s(1)/L))*u^2))/meff;
    thet*s(2)/Cp - s(3)/(R*Cp)];

%Call ode45 solver
[time, state_values] = ode45(sdot,tspan,IC);

%Extract individual values
x = state_values(:,1);
xdot = state_values(:,2);
y = state_values(:,3);

%plot x(t) and y(t)
figure(1)
clf
plot(time,x), title('Displacement'), xlabel('time(s)'), ylabel('displacement(m)')

figure(2)
clf
plot(time,y), title ('Voltage generation'), xlabel('time(s)'), ylabel('Output voltage(V)')

end

[Walter Roberson]

You need to call your function to get the output.

The output pretty much goes to infinity... but it does show up.

It is confusing that you folded your into other constants so we cannot be sure it is properly taken into account.

[Jayanth Palat]

Hello Walter!

Thank you for replying.

Even if i call the function the value only shows 0 for both x(t) and y(t) and then as you said, pretty much goes to infinity but, thats not the shape of the graph. The displacement should oscillate between positive and negative values.

What did you mean by foldin meff? That equation needs to be divided by meff, thats why I put it like that.

Also, how can I code the below part and put it in?

=∑.r=1,2…A=∑r=1,2….Ar (α)r

Thank you :)

[Walter Roberson]

Sorry, I did not notice the /meff before. I was expecting an implementation with the same terms in the same order as the equations.

sum(A .* (s(2)/U + 3*s(1)/L).^(1:length(A)))

[Jayanth Palat]

okay, I only need a1 and a3 but the output curve shape is not coming as expected. Is there any problem with with the extraction of the values from ode5?

a1 =2.3

a3 = -18

It seems as though the values extracted are only for the initial parts of tspan.

[Walter Roberson]

When I test, the values seem to be extracted without difficulty.

If you plot(time) you will see that it rises fairly steadily until roughly the last 30 points, and then it goes through the rest of the time from 22 to 60 within a relatively small number of points.

What is happening there is not it somehow omitting points. What is happening is that ode45() uses adaptive step sizes, and when it reaches that time, the slope of the curve has started to become steep enough that MATLAB feels comfortable taking large time steps.

[Jayanth Palat]

So how do I change the code to get the required output?

function [time, w, wdot, v] = bluff
meff = 0.0114; %Effective mass (kg)
ceff = 0.0212; %Effective damping (N/(m/s))
keff = 8.86;   %Effective stiffnes (N/m)
thet = 0.000019758; %Electromechanical coefficient (N/V)
Cp = 0.0000000157; %Capacitance
rho =1.225 ; %air density (kg/m3)
u = 10; %wind velocity(m/s)
stip = 0.0118; %exposed area of the bluff body(m2)
L = 0.15 ; %length of the beam
a1 = 2.3; a3 = -18 ; %empirical coefficient for the aerodynamic force calculation
R = 1050000; %load resistance (ohm)
q = R*Cp ;

%values
x0 = 0; xdot0 = 0; y0 = 0;
t0 = 0; tf = 60;

%Time span
tspan = [t0, tf];


%initial conditions
IC = [x0, xdot0,y0];

% [sdot] = g(t,s)
sdot = @(t, s) [s(2);
    -ceff*s(2)/meff - keff*s(1)/meff - thet*s(3)/meff + 0.5*rho*stip*u*u*(a1*((s(2)/u)+(1.5*s(1)/L))+a3*((s(2)/u)+(1.5*s(1)/L))^3)/meff ;
    thet*s(2)/Cp - s(3)/q];

%Call ode45 solver
[time, state_values] = ode45(sdot,tspan,IC);

%Extract individual values
w = state_values(:,1);
wdot = state_values(:,2);
v = state_values(:,3);

%plot x(t) and y(t)
figure(1)
clf
plot(time,w), title('Displacement'), xlabel('time(s)'), ylabel('displacement(m)')

figure(2)
clf
plot(time,v), title ('Voltage generation'), xlabel('time(s)'), ylabel('Output voltage(V)')

end

[Walter Roberson]

As a remote observer: we have no reason to know that the plot that is obtained is not in fact the correct plot for the equations.

[Walter Roberson]

You divide by Cp but Cp is about Pi / 2e8 . When you divide by that, values are magnified a lot.

You might have better success if y0 were a lot smaller.

[Jayanth Palat]

y0?

Walter, y0 = 0. Thats the initial condition.

Answer 1

[Lewis Fer]

Hello, I am having troubles solving a system of second order nonlinear equations with boundary conditions using MATALB

Here is the equations:

f''(t)=3*f(t)*g(t) -g(t)+5*t;

g''(t)=-4f(t)*g(t)+f(t)-7*t;

the boundary conditions are: f'(0)=0 et h'(o)=5;

g(0)=3 et h'(2)=h(2)

[Walter Roberson]

you should create your own Question for this purpose

[Lewis Fer]

I have post my question in this pub, but there is no answer

https://www.mathworks.com/matlabcentral/answers/853120-how-to-solve-a-system-of-second-order-nonlinear-differential-equations-with-boundary-conditions?s_tid=srchtitle#add_answer

[Lewis Fer]

https://www.mathworks.com/matlabcentral/answers/853120-how-to-solve-a-system-of-second-order-nonlinear-differential-equations-with-boundary-conditions?s_tid=srchtitle#add_answer