Integration drift with numerical simulations

조회 수: 11 (최근 30일)
Mishal Mohanlal
Mishal Mohanlal 2021년 10월 26일
댓글: Mathieu NOE 2021년 10월 27일
Hi Guys
I have solved a system of differential equations, whereby some of the variables are presenting integration drift in the results.
I have attached a spreadsheet with the results of one of the variables (theta) note that the time component was exluded to keep the file small. The time increments are 0.0001seconds with t starting at 0.
The plot of the results is presented below.
The inital conditions of the problem has theta = 0.1, which is correct when reviewing the data in the spreadsheet.
I have used the code below to process the data, however upon processing the data the inital condition is no longer true as the value is no longer 1. Note that the code does not present the data being read from the excel file as it already exists in a variable (theta).
Could someone please give me some advice on how one could accurately process the data?
T = detrend(theta,1);
figure(33)
plot(t,T)
Fs = 1/inc;
Fc = 1/(Fs/2);
[b a] = butter(4,Fc,'high');
T = filtfilt(b,a,T);
  댓글 수: 3
Mishal Mohanlal
Mishal Mohanlal 2021년 10월 26일
Hi Mathieu
Well let me rephrase, why does the value move off 0.1? I would have thought that the data at that point would stay at 0.1...
By doing the filtering I am basically messing up my data which is to be used for simulations
I would greatly appriciate some information so that I can understand.
Mathieu NOE
Mathieu NOE 2021년 10월 26일
hello Mishal
well doing filter can alter initial condition
if you use filter, this is only a forward time filter and you see the IC are kept the same so this is fine , but this way of doing filtering will cause a good amount of phase lag as will do any filter
if you use filtfilt , this applies the filtering both in forward and backward time axis, so this will cancel the phase lag and keep your signals time aligned , but the code cannot garantee that IC will be kept exact the same.
I will suggest you a solution in the answer section
see if it fills your needs and if you accept it

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

채택된 답변

Mathieu NOE
Mathieu NOE 2021년 10월 26일
here one possible solution
we want to have the signal not time distorted so we must keep with filtfilt even with the "bad" IC effects
a work around I suggest is to blend the unfiltered (raw) data and the filtered data within the first second, so that it will almost look like the raw data and then progressively match only the filtered data . So IC are kept , whatever the characteristics of your high pass filter
zoom on the start
overal behaviour
code :
theta = importdata('theta.txt');
theta = theta.data;
inc = 0.0001; % second
samples = length(theta);
t = inc*(0:samples-1);
Fs = 1/inc;
Fc = 1/(Fs/2);
[b, a] = butter(2,Fc,'high');
T = filtfilt(b,a,theta);
% blending the data on the first second
blend = ones(size(t))';
tmp = (0:inc:1);
blend(1:length(tmp)) = tmp;
T2 = blend.*T + (1-blend).*theta;
plot(t,theta,t,T,t,T2)
legend('raw','HP filtered','HP filtered and blended')
  댓글 수: 2
Mishal Mohanlal
Mishal Mohanlal 2021년 10월 27일
Hi Mathieu
Thanks for the assistance.
I am not sure if this is exactly what I am looking for but it does put me on the right track.
Mathieu NOE
Mathieu NOE 2021년 10월 27일
hello Mishal
ok - tx for accepting my submission
I have tried many other options, but none could really give me the right output
finally went to this trick as best available short term solution....

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

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Bartlett에 대해 자세히 알아보기

제품


릴리스

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by