Error function optimisation method to calculate the velocity of a travelling wave

Hi everyone, I'm trying to use an error optimisation method using interpolation and least quare method for a travelling wave solution. I have a mat file (SS) of size 10000*200 where time steps are defined by rows and node number is defined by columns. In this simulation, there are 10000 timesteps and 200 nodes. These values provide a travelling wave solution as follows.

To calculate the velocity of these waves, I'm trying to use the error optimisation method as given in the pseudo code below.

I tried coding that as follows, however, i'm having a trouble understanding and defining "c" in the anonymous function in the code.

v_t= zeros(Tmax-1);
sol_current = SS(1,:);
%x1= 1:Tmax/200:Tmax;
xq = linspace(0,1,200);
for n=1:length(Tmax)
    sol_future = SS(n+1,:);
    fun=@(c)...

end

It would be great if anyone can help me to code this method. Thank you!

Answer 1

[Torsten]

For each of the 200 times, extract the node number where your solution equals 0.5. Plot node number as a function of time. The slope of this (in your case most probably linear) curve gives the wave speed.

[UserCJ]

@Torsten This method didn't work as there were no exact node numbers for sol=0.5 for most of the time steps. So it returned ocsillations instead of an increasing graph.

[Torsten]

For each column, use interpolation to find the (not necessarily integer) position where your curve attains 0.5. The plot of position vs. time will be increasing throughout since the curves you posted don't intersect.

[UserCJ]

@Torsten could you please post the code chunk here?

[Torsten]

Assuming "x" is the 200x1 node array and "t" is the 10000x1 time array, you could do:

pos = zeros(10000,1);
for it = 1:10000
    vec = SS(it,:);
    pos(it) = interp1(vec,x,0.5);
end
plot(t,pos)

Maybe you could post a picture of your t-pos plot to look for the best way to get its slope.

[UserCJ]

Thanks @Torsten. I did it as follows; where N is the nodes array

pos = zeros(10000,1);
for it = 1:10000
    vec = SS(it,:);
    pos(it) = interp1(vec,N,0.5);
end
plot(t,pos)

But I'm getting the following error

Error using matlab.internal.math.interp1
Sample points must be unique.

Error in interp1 (line 188)
        VqLite = matlab.internal.math.interp1(X,V,method,method,Xqcol);

[Torsten]

In order to apply "interp1", the vector vec must be strictly monotonic. I hope that restricting the interpolation to values > 0.2 and < 0.8 will be successful to achieve this. Otherwise you have to include your data or experiment for yourself.

I assume N = 0:199, am I right ?

pos = zeros(10000,1);
for it = 1:10000
    vec = SS(it,:);
    idx = vec > 0.2 & vec < 0.8;    
    pos(it) = interp1(vec(idx),N(idx),0.5);
end
plot(t,pos)

[Star Strider]

In my approach to this problem (Answer, and a subsequent Comment with a slightly different implementation), I define a narrow range of indices to interpolate over (‘idxrng’ in that code) in each iteration of the loop. That usually works to eliminate the non-unique values problem, although I do not have the actual data mentioned here to experiment with.