Why is the resulting time-domain curve from Euler method and ode45 not stable?

조회 수: 1 (최근 30일)
I am writing a program to obtain the solution of an undamped free vibration system. At first, I used the Euler method and the resulting time-domain curve seemed normal (perfectly sinusoidal with the initial displacement as the amplitude). But when I checked the peaks using findpeaks, its value increased continuously (0.1, 0.1001, 0.1002, and so on). Then, I used the ode45 to do the same. The peaks of the time-domain curve decreased continuously from the initial displacement (0.1, 0.0999, 0.0998, and so on). Why is this happening? Am I doing something wrong?
  댓글 수: 2
Mischa Kim
Mischa Kim 2021년 1월 12일
Please attach your code so we can help.
Ni Made Ayu Sinta Dewi
Ni Made Ayu Sinta Dewi 2021년 1월 13일
Thank you for reading my problem. This is the code using Euler method.
m=1;
k=16;
kr=0.5;
ks=k*kr;
ke=(k*ks)/(k+ks);
c=0;
fs=100000;
dt=1/fs;
t=0:dt:20;
F1=0*sin(2*pi*t);
F=F1;
A=[0 1;-ke/m -c/m];
B=[0;1/m];
y=zeros(2,length(t));
y(1,1)=0.1;
y(2,1)=0;
for i=2:length(t)
dx=A*y(:,i-1)+B*F(:,i-1);
x=y(:,i-1)+dx*dt;
y(:,i)=x;
end
%FFT
L=numel(t)-1;
f=0:fs/L:fs/2;
Y=fft(y(1,:));
P2=abs(Y/L);
P1=P2(1:L/2+1);
P1(2:end-1)=2*P1(2:end-1);
subplot(2,1,1)
plot(t,y(1,:));
subplot(2,1,2)
plot(f,P1);
[pks_y,locs_t]=findpeaks(y(1,:));
The resulting peaks are shown here. As I have said previously, the amplitude is increasing continuously.
>> format long
>> pks_y
pks_y =
Columns 1 through 5
0.100007255462951 0.100014511451836 0.100021767966689 0.100029025007554 0.100036282574466
Columns 6 through 7
0.100043540667467 0.100050799286586
This is the code using ode45.
t0=0;
tf=20;
T=0.00001;
n=(tf-t0)/T;
seq=0:n;
tspan=T*seq;
x0=[0.1 0]';
[t,x]=ode45(@state_space,tspan,x0);
[ba,ko]=size(x);
N=ba;
T=tf/ba;
k=0:N-1;
f=k*(1/(N*T));
magF=2/ba*abs(fft(x));
subplot(211)
plot(t,x(:,1))
subplot(212)
N1=round(N/2000);
plot(f(1:N1),magF(1:N1,1),'r')
[pks,locs]=findpeaks(x(:,1));
function dx=odefun(t,x)
m=1;
k=16;
kr=0.5;
ks=k*kr;
ke=(k*ks)/(k+ks);
A=[0 1;-ke/m 0];
dx=A*x;
The resulting peaks are shown here. The amplitude is decreasing continuously.
>> format long
>> pks
pks =
0.100004887335947
0.099986382626947
0.099968120795277
0.099949910847763
0.099931713047562
0.099913520122880
0.099895330803077

댓글을 달려면 로그인하십시오.

채택된 답변

Mischa Kim
Mischa Kim 2021년 1월 13일
Hi Ni Made Ayu Sinta Dewi, there is nothing you are doing wrong. What you are seeing is expected behavior. This is because you are using a numerical solver (e.g. ode45) to solve the differential equation. Numerical solutions are not exact solutions, they are only approximations. However, you have some control regarding the accuracy of the solution. E.g. you can set tolerance levels:
options = odeset('RelTol',1e-13); % the smaller the tolerance, the more accurate the result
[t,x] = ode45(@state_space,tspan,x0,options);

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Acoustics, Noise and Vibration에 대해 자세히 알아보기

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by