# Find the velocity of a travelling wave like behaviour

5 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.
##### 2 CommentsShow 1 older commentHide 1 older comment
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
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!