Large error between exact analytic answer and numeric integration
1 次查看(过去 30 天)
显示 更早的评论
I am trying to solve the equation below in (1) and (2) using the numerical integrator (ode45, ode15s, etc.) within Matlab so I can compare it to the exact analytical steady-state solution from basic vibration analysis. This is a simplified problem to determine the cause of other issues I’m running into with other simulations.
(1) M= [3 1; 1 5],K=[7 -1; -1 10],F1= [1; 4],F2= [-1; 7]
(2 ) M(d^2x/dt^2)+Kx=F1 sin(αt)+F2cos(βt)
According to vibration theory, the exact steady-state analytical solution should be:
(3) x(t)=sin(αt)*(-Mα^2+K)^(-1)*F1+cos(βt)*(-Mβ^2+K)^(-1)*F2
Substituting in the matrices and performing the calculations, the result is:
(4) x(t)= [(-α^2-14)/(14α^4-67α^2+69); (-11α^2+29)/(14α^4-67α^2+69)] sin(αt)+[(12β^2-3)/(14α^4-67α^2+69); -22β^2+48)/(14α^4-67α^2+69))]cos(βt)
In my problem I assume α,β=1, so that the resulting solution is:
(5) x(t)=[0.8125; 1.125]sin(t)+[0.5625; 1.625]cos(t)
Approaching this problem using the numerical integrator (ode45, ode15s, etc.), the problem set up is as follows:
function dy = EoMExample(t,y)
% function to be integrated
%y(1) = x1 dy(1) = x1dot
%y(2) = x2 dy(2) = x2dot
%y(3) = x1dot dy(3) = x1dotdot
%y(4) = x2dot dy(4) = x2dotdot
%y(5) = x1dotdot dy(5) = x1dotdotdot
%y(6) = x2dotdot dy(6) = x2dotdotdot
%y(7) = x1dotdotdot
%y(8) = x2dotdotdot
dy = zeros(8,1);
dy(1) = y(3);
dy(2) = y(4);
dy(3) = (1/3)*(-y(6)-7*y(1)+y(2)+sin(t)-cos(t));
dy(4) = (1/5)* (-y(5)+y(1)-10*y(2)+4*sin(t)+7*cos(t));
dy(5) = y(7);
dy(6) = y(8);
Where the driver is:
%% Driver for integrating equations of motion
runtime = 30; %[sec] Total desired runtime for each
free stream velocity
%Initial Conditions
x10 = 0;
x20 = 0;
x1dot0 = 0;
x2dot0 = 0;
x1dotdot0 = 0;
x2dotdot0 = 0;
x1dotdotdot0 = 0;
x2dotdotdot0 = 0;
% Integrate EoM for the specified timespan
tspan1 = [0 runtime];
x_init = [x10; x20; x1dot0; x2dot0; x1dotdot0; x2dotdot0; x1dotdotdot0; x2dotdotdot0];
options = odeset('abstol',1e-8,'reltol',1e-5,'InitialStep',1e-4,'MaxStep',1e-3);
[t, x] = ode45(@(t,y) EoMExample(t,y), tspan1, x_init, options);
If I run this driver and plot the results for x1 and x2, it is very chaotic and does not reach a steady-state value. The two states are not simple sinusoids as represented by the analytic solution (eqn. (5)). My question is how there could be such a major discrepancy? I have tried running the integrator for large values to see if steady-state is achieved, but it never happens. I’ve tried stiff and non-stiff integrators with the same result from both. Is there some fundamental aspect of this problem I have overlooked? Are the numerical integrators within Matlab simply not able to handle coupled 2nd order differential equations? Any and all feedback is much appreciated.
0 个评论
回答(2 个)
Mark Shore
2012-2-24
I rarely work with differential equations but you should check your basic equation again.
For some reason you have a 1D harmonic oscillator with NO damping term and two sinusoidal driving forces. Normally there should be a dx/dt term (even a small one) in your equation.
Without a damping term there is nothing to stop the amplitudes from growing, generating numerical instability and never reaching a steady state.
Mike Hosea
2012-2-28
Your formulation has many problems. I don't see that you're using a mass matrix with ODE45, and you also have not inverted M, at least not correctly. Why 8 components? You're setting the 2nd and 3rd derivatives to zero. The steady state solution does not have that property, so obviously this formulation is incorrect. When I substitute a correct formulation of the problem (with only 4 components) and begin the integration on the steady-state solution, it stays there.
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!