I will try to show how we can solve an Ordinary Differential Equation using Matlab solver ODE45. However, with this problem i will use GNU Octave environment to solve. First, i will model a car suspension damping system equations then solve it.
In my previous post analogies-between-elements-in-electrica ... stems-5072, I modeled Mass, Spring, Damper System with its Electric circuit equivalent. So, in this figure a mechanical damper is presented and i will model and solve it at one initial condition as well as different initial conditions for displacements and velocity.
Balance of forces in Mass, spring, damper system:
ODE45 format:
- [t,x] = ode45(@FunctionName, tspan, [xInitial; xdotInitial])
- FunctionName = [x(2), RightHandSideofAFunction]
Where
m-mass
D-Damping Ratio
K-Spring constant
F-Force
Let
then and
Replace the balance of Force equation with these new parameters and make the subject;
This plot shows that the Suspension system overshoots [Displacement and Velocity ] initially and becomes stable after some time on different initial conditions [from -1 to 1] , I think with the Phase Plane we might have a limit circle if i extend the initial conditions for displacement. This system has local equilibrium point as seen in the Phase Plane plot.
Asymptotically Stable System
What determines the quality of the systems, is it mass, Force, Damping Coefficient, or Spring Constant. I will modify the code to plot for different Damping coefficient and i will see the differences.
Global Stable system
So, in choosing the type of material to make a suspension system, the parameters above are taken into account to have a durable "shokapu"
@Admin please help me to put the codes below in a good listing, its .m file
- function kkoo = damper(t,x)
- m = 1; %mass
- F = 10; %Force
- D = 0.01; %Damping Coefficient
- K = 10; %Spring Constant
- tzr = x(2); %let xdot =x(2)
- suka = 1/m.*(F - D*x(2)-K*x(1)); % Right hand side of the ODE
- kkoo = [tzr; suka];
- endfunction
- function venom
- t = 0:0.1:20;
- initial_x = 0; % different initial condition, you can put any range say [-1:1]
- initial_xdot = 0;
- for j = 1:length(initial_x)
- for k = 1:length(initial_xdot)
- cond3 = initial_xdot(k);
- cond2 = initial_x(j);
- [t,x] = ode45(@damper, t , [cond2 cond3] );
- figure 1
- h(1) = subplot(2,2,1);
- h(2) = subplot(2,2,2) ;
- h(3) = subplot(2,2,[3 4]);
- set(h(1),'NextPlot','add');
- set(h(2),'NextPlot','add');
- set(h(3),'NextPlot','add');
- grid on
- plot(h(1),t,x(:,1));
- xlabel('Time t');
- ylabel('Displacement x');
- title('Car suspension model displacement');
- grid on
- plot(h(2),t,x(:,2));
- xlabel('Time t');
- ylabel('Velocity-xdot');
- title('Car suspension model Velocity');
- grid on
- plot(h(3),x(:,1),x(:,2));
- xlabel('Displacement x');
- ylabel('Velocity-xdot');
- title('Car suspension model Phase Plane');
- endfor
- endfor
- endfunction