Problem in finding peak, how to find the first peak?

Hi, I'm trying to plot the P value (corresponding to a peak in a P-umax plot) in function of "s", a variable linspaced between 0.1 and 2.

the code that produces this plot P-umax for "s=1" is the following:

clc;
clear;
close all;
R=20;
L=3*75;
s=1;
ni=0.3;
E=210000;
w=2;
K=0.5;
%

E1=E/(1-(ni^2));
I=(2*pi*R*(s^3))/12;
k=2*pi*((E*s)/R);
P = linspace(400000,9000000);
y_max_values = zeros(size(P));
for i=1:length(P)
    f = @(z,u) [u(2); u(3); u(4); (-1/(E1*I))*((P(i)*u(3))+(k*u(1)))];
    bc = @(ua, ub) [ua(1); ua(2)-0.00000001; ub(1); ub(2)+0.00000001];
    zmesh = linspace(0, L, 100000);
    iguess= @(z) [0; 0; 0; 0];
    solinit = bvpinit(zmesh, iguess);
    Sol= bvp4c(f, bc, solinit);
    
    x = Sol.x;
    y = Sol.y;
    idx = (x > 50) & (x < 180);
    x_new = x(idx);
    y_new = y(1, idx);
    
    y_max_values(i) = max(y_new);
end
figure
hold on
grid on
plot(y_max_values, P)
legend ('cilindro 1 w0')
xlabel ('u_max [mm]')
ylabel ('P [N]')

So I need to find in the plot, the P value corresponding to the first peak of the following P-umax plot, where the umax value moves from zero in the first time. (8x10^5 in this case), and repeat this procedure for every "s".

The P value corresponding to peak I'm trying to find has to be assigned to "Pcr"

I've tried to write something like this but I don't know what to put in the Pcr findpeak function

clc;
clear;
clear;
%------------------------------------------------------------------------
ni=0.3;
L=50;
R=20;
E=210000;
E1=E/(1-ni^2);
%-----------------------------------------------------------------------

s= linspace(0.1,2,10);
I = zeros(size(s));
k = zeros(size(s));
Pcr = zeros(size(s));
%%
P=linspace(7800,3300000,10);
y_max_values=  zeros(size(P));
for i=1:length(s)
    for j=1:length(P)
     f = @(z,u) [u(2); u(3); u(4); (-1/(E1*I(i)))*((P(j)*u(3))+(k(i)*u(1)))];
     bc = @(ua, ub) [ua(1); ua(2)-0.00000001; ub(1); ub(2)+0.00000001];
    zmesh = linspace(0, L, 100000);
    iguess= @(z) [0; 0; 0; 0];
    solinit = bvpinit(zmesh, iguess);
    Sol= bvp4c(f, bc, solinit);
    
    x = Sol.x;
    y = Sol.y;
    idx = (x > 50) & (x < 180);
    x_new = x(idx);
    y_new = y(1, idx);
    
     y_max_values(j) = max(y_new);
    Pcr=findpeak (?) %what I should put here?
    end
    
end
%%
figure
hold on
grid on

plot(s, Pcr)

xlabel ('s [mm]')
ylabel ('Pcr [N]')

Could someone help me out ?

Is the the code structure right for the aim I need to?

[Ameer Hamza]

Can you mark the peak in the image you are trying to find? Also, the code takes a long time to run; it will be easier to suggest a solution if you attack vector y_max_values and P in a .mat file?

[federico midei]

Thank you for your interest! Really appreciate it

here it is:

how should I do the attacking vector thing?

[federico midei]

So...do you have any suggestion please? I can’t attach P and y in a file ‘cause they’re part of the computation themselves, I can’t do a file.m for them

[Star Strider]

I tried to run your code (in R2020b) and could not, so I gave up.

The error was:

Error using bvp4c (line 251)
Unable to solve the collocation equations -- a singular Jacobian encountered.

.

[federico midei]

I solved the "singular" jacobian problem, but now I get that:

Unable to perform assignment because the indices on the left side are not compatible with the size of the right side.

Error in plot_s_Pcr (line 50)
     y_max_values(j) = max(y_new);
     

I don't understand why, because in the code running without the "for J" cycle I don't get this error.

here is my code now

     clc;
clear;
clear;
%------------------------------------------------------------------------
ni=0.3;%                            -coeff poisson
L=50;%mm                            -lunghezza tubo
R=20;%mm                            -raggio medio
E=210000;%MPa %acciaio              -modulo young
E1=E/(1-ni^2);%                     -modulo young dilataz laterale impedita
I=(2*pi*R*(s.^3))/12;%       
k=2*pi*((E.*s)/R);%       



P=linspace(7800,3300000);
s= linspace(0.1,2,10);
for i=1:length(s)
    for j=1:length(P)
        %Ii=I(i);
        %ki=k(i);
        %Pj=P(j);
     f = @(z,u) [u(2); u(3); u(4); (-1/(E1*I(i)))*((P(j)*u(3))+(k(i)*u(1)))];
     bc = @(ua, ub) [ua(1); ua(2)-0.00000001; ub(1); ub(2)+0.00000001];
    zmesh = linspace(0, L, 100000);
    iguess= @(z) [0; 0; 0; 0];
    solinit = bvpinit(zmesh, iguess);
    Sol= bvp4c(f, bc, solinit); %sol for diff equation
    
     % this takes the max y value in sol.y and stores it in y_max_values
     x = Sol.x;
     y = Sol.y;
    idx = (x > 50) & (x < 180);
    x_new = x(idx);
    y_new = y(1, idx);
    
     y_max_values(j) = max(y_new);
    %    I NEED TO FIND THE FIRST PEAK IN EVERY P-y_max_values PLOT
    %Pcr=findpeaks();
    end
    
end

figure
hold on
grid on
plot(s, Pcr,'o')
legend ('attended by calc.')
xlabel ('s [mm]')
ylabel ('Pcr [N]')

Is someone able to debug this please?

[Star Strider]

The easiest way to solve that problem is to save the ‘max(y_new)’ (that is likely a row vector) is either:

y_max_values(j,:) = max(y_new);

or as a cell array, if the size is changing:

y_max_values{j} = max(y_new);

.

[federico midei]

Thank you!

found the problem, my bad. It was an error on the L vaue that was <50 so the idx couldn't work.

Do you have any advice on finding the peaks?

as I wrote in the code : %I NEED TO FIND THE FIRST PEAK IN EVERY P-y_max_values PLOT

%Pcr=findpeaks();

[Star Strider]

My pleasure!

Please edit the code. When I try to run the version you posted, I am still getting the error, and it is taking a while to run.

EDIT — I stopped it after it ran for several minutes with no results. Please post a working version of your code.

[federico midei]

ok sorry.

Here a recap on what I need to do:

so In the loop "j" I have the y_max_values that varies in function of P. If the "s" variable was a constant value i could do a plot P(on the x axys) and y_max_values(on the y axys).

Now, for every "s" value (so in every "i" loop) I need to take the P value corresponding to the first value of y_max_values that results >0, store the value in "Pcr" so called.

finally I have to plot every Pcr value in function of every s value.

Here you have the code:

clc;
clear;
close all;

R=20;%mm                                        -raggio medio
L=3*75;%mm                                      -lunghezza tubo

ni=0.3;%                                        -coeff poisson
E=210000;%MPa                                   -modulo young
w=2;%MPa                                        -pressione int/ext
K=0.5;%                                         -incasto incastro


%equazioni derivanti dall'equilibrio radiale:

E1=E/(1-(ni^2));%                     -modulo young dilataz laterale impedita



s= linspace(0.1,2);

I=(2*pi*R*(s.^3))/12;%                 -momento di inerzia cilindro
k=2*pi*((E.*s)/R);%                    -rigidezza molla fondazione
Pcr = zeros(size(s));

P=linspace(400000,3300000);
y_max_values= zeros(size(P));
for i=1:length(s)
    for j=1:length(P)

        f = @(z,u) [u(2); u(3); u(4); (-1/(E1*I(i)))*((P(j)*u(3))+(k(i)*u(1)))];
        bc = @(ua, ub) [ua(1); ua(2)-0.00000001; ub(1); ub(2)+0.00000001];
        zmesh = linspace(0, L, 100000);
        iguess= @(z) [0; 0; 0; 0];
        solinit = bvpinit(zmesh, iguess);
        Sol= bvp4c(f, bc, solinit);
    
  

    
     x = Sol.x;
     y = Sol.y;
    idx = (x > 50) & (x < 180);
    x_new = x(idx);
    y_new = y(1, idx);
    %y_max_values{j} = max(y_new);
    %y_max_values(j,:) = max(y_new);
    y_max_values(j) = max(y_new);
    
    
    %here I need to do the peaks thing that leads to define the Pcr variable
    
    end
end

figure
hold on
grid on
plot(s, Pcr,'o')
legend ('attended by calc.')
xlabel ('s [mm]')
ylabel ('Pcr [N]')

[Star Strider]

The problem is that your code takes forever to run on my desktop machine (relatively new AMD Ryzen 9 3900 12-Core CPU). I simply cannot tie it up for several hours.

Put this line just after the loops (forget about findpeaks for the time being), then run your code and let it finish:

save('federico midei_20201008.mat', 'y_max_values','P','s')

Add or delete from this line whatever variables you want (I included the ones that appear to be most important). Then run your code and attach the federico midei_20201008.mat file to your original post or to a Comment.

We can then explore those results.

[federico midei]

I’ ve stopped it ‘cause got it run for 8+ hours...in these case there is something to do in order to speed up the code?

these are the values for vectors s and p with length 10

Answer 1

[Star Strider]

I do not see any peaks. It is just sort of L-shaped lines.

That aside, if you do have data with peaks and you only want the first one (with the lowest value of ‘P’), and since your data appears to have no noise:

[pks,locs,wdth,prom] = findpeaks(P, 'NPeaks',1)

will llikely do what you want. It will return only the data for the first peak, which should be the one at about in the image you posted. (I experimented with this with the sinc function that definitely has peaks.)

To be certain you have the one you want, plot it as:

figure
plot(u_max, P)
hold on
plot(u_max(locs), pks, '+r')
hold off

If you want the one on the top of the plot instead (at about ), try this:

[pks,locs,wdth,prom] = findpeaks(P)

figure
plot(u_max, P)
hold on
plot(u_max(locs(end)), pks(end), '+r')
hold off

You can then address them the same way to get their values and locations.