Find the velocity of a travelling wave like behaviour

4 次查看（过去 30 天）
UserCJ
UserCJ 2023-7-31
评论： Bruno Luong 2023-9-11
Hi everyone, I have a travelling wave solution drawn at each time step for a set of data obtained from a previous simulation.
Tmax = 10000;
m =load('solution.mat');
SS = m.sol;
N = m.N1;
for i = 1:50:Tmax
Xval = SS(i,N);
hold on
plot(N,Xval,'Linewidth',2);
end
hold off
and below is the out put I got from the above code where x axis represents a distance (N) from 0 to 200 and y axis takes values between 0-1. .Each coloured line is drawn at a different time step.
Next, I'd like to calculate the velocity at each time step when y=0.5. I.e. I want to identify the x value when y=0.5 for each colour line and calculate the velocity as . The main question here is sometimes there are no exact points in Xval with value 0.5 or no exact N value associated to 0.5.
I'd be happy if someone could help me in creating a loop for this calculation.
Thanks in advance!
  2 个评论
Image Analyst
Image Analyst 2023-7-31
编辑：Image Analyst 2023-7-31
Ran in:
Tmax = 10000;
m =load('solution.mat');
Error using load
Unable to find file or directory 'solution.mat'.
SS = m.sol;
N = m.N1;
for i = 1:50:Tmax
Xval = SS(i,N);
hold on
plot(N,Xval,'Linewidth',2);
end
hold off
Where is it? You forgot to attach it!
What is velocity? What variable? What are the units of the x axis and y axis?
Make it easy for us to help you, not hard.
UserCJ
UserCJ 2023-8-1
@Image Analyst It is a large file, so it doesn't let me upload it even after compressing.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Bruno Luong
Bruno Luong 2023-8-1
Ran in:
%s =load('solution.mat');
%SS = s.sol;
% Create fake N and SS
N = -20:0.1:20;
v = arrayfun(@(x) x/(sqrt(1+x^2)), (1:40)/10)
v = 1×40
0.0995 0.1961 0.2873 0.3714 0.4472 0.5145 0.5735 0.6247 0.6690 0.7071 0.7399 0.7682 0.7926 0.8137 0.8321 0.8480 0.8619 0.8742 0.8849 0.8944 0.9029 0.9104 0.9171 0.9231 0.9285 0.9333 0.9377 0.9417 0.9454 0.9487
x0 = -17 + cumsum(v);
[X,X0] = meshgrid(N,x0);
SS = 1-(tanh(X-X0)+1)/2;
% Determine cross point
yx = 0.5;
ys = SS'-yx;
[m,n] = size(ys);
[i,j] = find(ys(1:end-1,:).*ys(2:end,:) <= 0);
% linear interpolation
i1 = i + (j-1)*m; ss1 = ys(i1);
i2 = i1 + 1; ss2 = ys(i2);
w = ss2./(ss2-ss1);
x = N(:);
ix = nan(1,n);
ix(j) = w.*x(i) + (1-w).*x(i+1);
vest = gradient(ix)
vest = 1×40
0.1961 0.2417 0.3294 0.4093 0.4809 0.5440 0.5991 0.6468 0.6880 0.7236 0.7541 0.7804 0.8032 0.8229 0.8400 0.8550 0.8681 0.8795 0.8897 0.8986 0.9066 0.9137 0.9201 0.9257 0.9309 0.9356 0.9397 0.9436 0.9470 0.9502
figure
plot(N, SS.');
hold on
plot(ix, yx+zeros(size(ix)), '+r', 'Markersize', 10)
figure
plot(vest)
hold on
plot(v)
legend('v estimate', 'v', 'location', 'best')
  7 个评论
显示 5更早的评论
Bruno Luong
Bruno Luong 2023-9-11
编辑：Bruno Luong 2023-9-11
Not sure why you need to filter out the oscillation, if the physical model and simulation is correct, then the oscillation is the physics reality.
If not then better figure out why the oscillation artefact occurs.
But here I do 2 fits with polynomial and spline. Spline has better behavior according to my eyes. I also fits on peaks+valleys only, since your data look skewed (the peaks is sharper) so that the fit goes somewhat close to the middle of alternating peak and valeys.
load('matlab.mat')
n = length(vest);
x = 1:n;
for fitfunction = {'poly' 'spline'}
b = islocalmin(vest) | islocalmax(vest);
switch fitfunction{1}
case 'poly'
[P, ~, mu] = polyfit(x(b), vest(b), 4);
vestsmooth_p = polyval(P,(x-mu(1))/mu(2));
case 'spline'
% Need file here https://www.mathworks.com/matlabcentral/fileexchange/25872-free-knot-spline-approximation
pp = BSFK(x(b), vest(b));
vestsmooth_s = ppval(pp,x);
end
end
figure
h1=plot(x,vest);
hold on
h2=plot(x,vestsmooth_p,'k','LineWidth',1);
h3=plot(x,vestsmooth_s,'r','LineWidth',1);
legend('data','polynomial','spline')
Bruno Luong
Bruno Luong 2023-9-11
Ran in:
For completness this is matlab spline fitting. It oscillates more than I prefer.
load('matlab.mat')
n = length(vest);
x = 1:n;
for fitfunction = {'poly' 'spline'}
b = islocalmin(vest) | islocalmax(vest);
xb = reshape(x(b),[],1);
y = reshape(vest(b),[],1);
switch fitfunction{1}
case 'poly'
[P, ~, mu] = polyfit(xb,y, 6);
vestsmooth_p = polyval(P,(x-mu(1))/mu(2));
case 'BSFK'
% Need file here https://www.mathworks.com/matlabcentral/fileexchange/25872-free-knot-spline-approximation
pp = BSFK(xb, y);
vestsmooth_b = ppval(pp,x);
case 'spline'
N = 12;
[xmin, xmax] = bounds(xb);
xi = linspace(xmin, xmax, N);
B = zeros(numel(xb),N);
for k=1:N
yi = accumarray(k,1,[N,1]);
B(:,k) = spline(xi, yi, x(b));
end
c = B \ y;
vestsmooth_s = spline(xi, c, x);
end
end
figure
h1=plot(x,vest,'g');
hold on
h2=plot(x,vestsmooth_p,'k','LineWidth',1);
%h3=plot(x,vestsmooth_b,'r','LineWidth',1);
h4=plot(x,vestsmooth_s,'r','LineWidth',1);
legend('data','polynomial','spline')

请先登录，再进行评论。

更多回答（1 个）

KSSV
KSSV 2023-8-1
I will suggest to use the function InterX.
N = 1000 ;
x = linspace(0,200,N) ;
y = rand(size(x)) ;
L1 = [x;y] ;
L2 = [x; repelem(0.5,1,N)] ;
P = InterX(L1,L2) ;
You got the points you want in P. You can do what you want
.

类别

在 Help Center 和 File Exchange 中查找有关 Polynomials 的更多信息

标签

产品


版本

R2022b

Community Treasure Hunt

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

Start Hunting!

Translated by