# Find the velocity of a travelling wave like behaviour

7 views (last 30 days)
UserCJ on 31 Jul 2023
Commented: Bruno Luong on 11 Sep 2023
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;
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.
Image Analyst on 31 Jul 2023
Edited: Image Analyst on 31 Jul 2023
Tmax = 10000;
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?
UserCJ on 1 Aug 2023
@Image Analyst It is a large file, so it doesn't let me upload it even after compressing.

Bruno Luong on 1 Aug 2023
%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 = 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')
Bruno Luong on 11 Sep 2023
Edited: Bruno Luong on 11 Sep 2023
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.
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 on 11 Sep 2023
For completness this is matlab spline fitting. It oscillates more than I prefer.
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')

KSSV on 1 Aug 2023
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
.

### Categories

Find more on Shifting and Sorting Matrices in Help Center and File Exchange

R2022b

### Community Treasure Hunt

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

Start Hunting!