I can't figure out weither the problem resides in your model or in the computations, but when I had (many years ago) to use discrete convolution to approximate continuous convolution, I noticed that the best way to obtain coherent results is to use the composed midpoint rule to approximate the integral
h(t)=\int_0^t f(\tau) g(t-\tau) d\tau. For two causal signals t->f(t)], t-> g(t) with respective support [0,N*δ] and [0,M*δ], sampled in the middle of each [i*δ,(i+1)*δ] interval, define the discrete sequences f_i = f((i-1/2)*δ), i=1...N g_i = g((i-1/2)*δ), i=1..M then define h_i = h((i*δ) for i=0...N+M by h_0 = 0 h_i = δ*\sum_{k=1}^i \tilde f_{i-k+1}*g_k, i = 1...N+M-1 h_{N+M} = 0 h_i for i = 1...N+M-1 is the approximation of the integral by the composed midpoint rule and is directly given in Scilab by δ*conv(f,g) where f and g are defined above. Here is an example: function y=f(t) y=sin(2*%pi*t); endfunction function y=g(t) y=sin(%pi*t); endfunction δ=0.01; N=50; // t->f(t) has support [0,0.5] M=100; // t->g(t) has support [0,1] h=[0 δ*conv(f((1/2:N)*δ), g((1/2:M)*δ)) 0]; th=(0:N+M)*δ; // t->h(t) has support [0,1.5] tf=(0:N)*δ; tg=(0:M)*δ; clf plot(tf,f(tf),tg,g(tg),th,h) legend f g f*g S. Le 03/02/2023 à 11:52, Heinz Nabielek a écrit : This is my latest code version: the 'convoluted secondary failure fraction' is nicely below primary failure, but seems too low. Heinz m=2; // Weibull modulus in mechanism #1 k=1E-7; // corrosion rate(s-1) in mechanism #1 n=500; // number of time steps in hourly intervals t=(1:n)'; // timesteps in hours ts= 3600*t; // timesteps in seconds PHI= 1 - exp(-((k*ts)^m)); // failure fraction in mechanism #1 deriv=m*k*((k*ts)^(m-1)) .* exp(-((k*ts)^m)); //derivative of PHI phi=3600*deriv; //derivative of PHI integrated over an hour f= 1 - exp(-ts/3E5); // "pure" failure fraction in mechanism #2 F=(convol(phi,f)/n)(1:n); //"convoluted" failure fraction mechanism #2 plot("nl",t, [PHI phi f F'],'-');xgrid; legend('PHI primary failure','phi primary failure rate','pure secondary failure fraction','convoluted secondary failure fraction',4); xlabel('time (hours)'); ylabel('failure fraction/ failure rate'); On 03.02.2023, at 11:24, Heinz Nabielek <heinznabie...@gmail.com><mailto:heinznabie...@gmail.com> wrote: On 03.02.2023, at 11:13, Stéphane Mottelet <stephane.motte...@utc.fr><mailto:stephane.motte...@utc.fr> wrote: Thanks for the code. Just a remark on the notations, you should write : F(T) = Int_{0}^{T} PHI(t) . f(T-t) . dt i.e. not F(t) since t is mute. However, you should pay attention to the delay notion associated with convolution and the relationships between discrete convolution and continuous convolution. I am not sure that the output of conv matches with a given discretization of the integral above. Maybe rectangle method, but I am not sure at all. Anyway, you should have F(0)=0 which does not seem to be the case in your graph. F(t), of course. But the formula is fundamentally wrong and one should have seen it already from a dimensional consideration. Correct should be F(t) = Int_{0}^{T} (d PHI(t)/dt) . f(T-t) . dt Hope you can help between conv and convol and possible something else. Heinz Dear SciLab Friends: I have an object consisting of many (~10,000) small components that can fail in a statistical way during long-term operation at extreme conditions. My primary failure model is described by PHI(t) going monotonically from zero to one at times from t=0 to T. In the computer, this is realized in n timesteps. Mechanism 2: There is a secondary failure mechanism that starts only when a component has previously failed by mechanism 1 and occurs after some delay. The delay function is f(t) and final expression for the evolution of the mechanism 2 failure is given by: F(t) = Int_{0}^{T} PHI(t) . f(T-t) . dt a classical convolution Faltungsintegral. In Scilab, I write F= conv(PHI, f, "same")/ n; and the resulting function F looks reasonable for most of the time. Under some conditions, however, I can have F>PHI at early stages, but this is physically impossible. Because it comes later, F must always be below PHI. What do I do wrong? Heinz _________________________ PS: massively simplified code below...... k=1E-7; // corrosion rate(s-1) in mechanism #1 n=100; // number of timesteps t=(1:n)'; // time in days ts= 3600*t; // time in seconds PHI= 1 - exp(-((k*ts)^2)); // failure fraction in mechanism #1 plot("nl", t, PHI, 'r--'); f= 1 - exp(-ts/3E5); // "pure" failure fraction in mechanism #2 F=conv(PHI,f,"same")/length(t); //"convoluted" failure fraction mechanism #2 plot(t,F,'b--'); legend('failure fraction #1','subsequent failure fraction #2',4); xlabel('time (hours)'); ylabel('failure fraction'); title('Component failure evolution #1 is followed by #2'); -- Stéphane Mottelet Ingénieur de recherche EA 4297 Transformations Intégrées de la Matière Renouvelable Département Génie des Procédés Industriels Sorbonne Universités - Université de Technologie de Compiègne CS 60319, 60203 Compiègne cedex Tel : +33(0)344234688 http://www.utc.fr/~mottelet This email and any attachments are intended solely for the use of the individual or entity to whom it is addressed and may be confidential and/or privileged. If you are not one of the named recipients or have received this email in error, (i) you should not read, disclose, or copy it, (ii) please notify sender of your receipt by reply email and delete this email and all attachments, (iii) Dassault Systèmes does not accept or assume any liability or responsibility for any use of or reliance on this email. Please be informed that your personal data are processed according to our data privacy policy as described on our website. Should you have any questions related to personal data protection, please contact 3DS Data Protection Officer https://www.3ds.com/privacy-policy/contact/
_______________________________________________ users mailing list - users@lists.scilab.org Click here to unsubscribe: <mailto:users-unsubscr...@lists.scilab.org> https://lists.scilab.org/mailman/listinfo/users