Find the velocity of a travelling wave like behaviour

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!

[Image Analyst]

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]

@Image Analyst It is a large file, so it doesn't let me upload it even after compressing.

Answer 1

[Bruno Luong]

%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')

[UserCJ]

What does

[i,j] = find(ys(1:end-1,:).*ys(2:end,:) <= 0);

in cross point detection mean? Could you please explain the linear interpolation used in here? Thanks!

[Bruno Luong]

ys(1:end-1,:).*ys(2:end,:) <= 0

is true when ys changes the sign (or one of them is 0) meaning when ys cross 0, or in other word SS cross yx := 0.5.

The linear interpolation computes where the straight line segment L that links the point P=(N(i),ys(i)) to the next point Q=(N(i+1),ys(i+1) cross the horozontal axis y=0, using barycentric parametrization of the line.

L(w1) := w1*P +w2*Q, with w1+w2 = 1;

If you rename w = w1 and solve for w the equation y(w) = 0

0 = w*ys(i) + (1-w)*ys(i+1)

you'll find

w = ys(i+1) / (ys(i+1)-ys(i))

this is the implementation in the code

w = ss2./(ss2-ss1);

The x-coordinates of the crossing is then

w.*N(i) + (1-w).*N(i+1)

This kind of formula probably written in a tone of books, but simple enough of derive every time one needs it.

[UserCJ]

@Bruno Luong This is what I got from your method. I just wanted to eliminate oscillations in this plot.

[Bruno Luong]

@UserCJ Please save vest (the one you plot) in a mat file and upload as attached file here

[UserCJ]

@Bruno Luong see attached!

[Bruno Luong]

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]

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')

Answer 2

[KSSV]

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

.