How to obtain Duhamel's integral?

USERPC01 · Регистрация: 29.06.2018

Студворк — интернет-сервис помощи студентам

My file (C++):

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
 #include <stdio.h>
 
#include <math.h>
#include <conio.h>
 
 
double p(double tau)
{
/*dU/dt, t=tau;*/
return -20000*exp(-200*tau);
}
double R=1000;
double C=0.00000001;
double k=0.25;
 
 
double h(double t ) {
 
 
/*
double R=1000;
double C=0.00000001;
double k=0.25;
double T=R*C;
 
return k*(1-exp(-t/T);
*/
 
return  0.005-0.0025*exp(-100*t);
}
 
double f(double tau, double t  )
{
 
 
return p(tau)*h(t-tau);
 
}
 
 
double trap(double a,double b,double eps, double(*f) (double , double ))
{
       
double h1,s,s0,s1,sn, t;
int i,n;
s=1;
sn=101;
n=4;
s0=(f(a,t)+f(b,t))/2;
s1=f(((a+b)/2),t);
while (fabs(s-sn)>eps){
sn=s;
h1=(b-a)/n;
for (i=0; i<n/2; i++)
s1+=f((a+(2*i+1)*h1),t);
s=h1*(s0+s1);
n*=2;
}
return s;
}
 
double DuhamelInt( double t)
{
double a,b,er,eps;
double f(double, double),s,trap(double,double,double,double(*) (double,double)) ; 
eps=0.0000001;
 
return trap(0,t,eps,f);
}
 
void main()
{
FILE *fp;     
double  time,Tend,step,s;
printf ("\n s(t)= trap(0,t,eps,f); \n f(tau)=p(tau)*h(t-tau) ");
printf ("\n Input Tend, ms (120) ");
scanf("%lf",&Tend);
Tend=Tend*0.001;
printf ("\n Input step ,ms,(0.01) ");
scanf("%lf",&step);
 
 step=step*0.001;
 
 
if ((fp=fopen("Duhamel.txt","a"))==NULL)  printf ("\n File opening error");
 
fprintf (fp,"\n s(t)= trap(0,t,eps,f); \n f(tau)=p(tau)*h(t-tau) ");
fprintf (fp,"\n  Tbeg=0 , Tend=%le s ",Tend);
fprintf (fp,"\n   step=%le s ", step);
fprintf (fp,"\n   p(tau)=-20000*exp(-200*tau); "); 
fprintf (fp,"\n   h(t)=0.005-0.0025*exp(-100*t); \n");
 
for(time=0; time<=Tend; time=time+step)
{  
s=DuhamelInt(time);
printf("\n t= %le s u(t)=%le ",time,s) ;
fprintf (fp,"\n t= %le s u(t)=%le ",time,s);
}
fclose(fp);
getch() ;
}

How to obtain Duhamel's integral using Delphi?

Добавлено через 20 часов 42 минуты
Как использовать механизм задания функции , аналогичный

C++
1
 double(*f) (double , double ))

в строке

C++
1
 double f(double, double),s,trap(double,double,double,double(*) (double,double))

для гибкости , модульности, многоразовости и универсальности программы, чтобы не писать заново под каждую функцию ?

Добавлено через 17 минут
Еще один дополнительный вопрос : можно ли работать с парсером trap как с библиотекой trap.dll, многократно обращаясь к нему как к подпрограмме на С++?
Как передавать в него аргументы , функции на случай вычисления двумерного массива из этих функций , каждая со своими параметрами в рамках некоторого класса в зависимости от массива функций -моделей (или их параметров, например массив моделей h[i][j](t), f[i][j](tau[i][j],t) ,p[i][j](tau[i][j]) ), с процедурой-методом (Public ) их редактирования как метод класса? Как использовать поинтеры и динамические массивы, как лучше к ним приделать деструктор при уничтожении класса массива редактируемых адаптивных парсеров импульсных характеристик и передаточных функций (y[i][j](t)= f[i][j](x[i][j](t),t, params_array[i][j],functions[i][j] ))?

Добавлено через 13 минут
( в смысле функций вычисления отклика (y[i][j](t)= f[i][j](x[i][j](t),t, params_array[i][j],functions[i][j] )) и адаптивных параметрически зависимых моделей функций импульсных характеристик )?

@Kitayets · 06.07.2018, 01:58

функция DuhamelInt выглядит странной. похоже на предварительное описание функций (типа forward в delphi), что имело бы смысл, если бы функции f и trap были бы описаны ниже или вообще бы линковались из другого объектного модуля, в данном случае строки 64 (кроме переменной eps) и 65 можно вообще выкинуть. зачем объявлять переменные которые не используются??? явно код откуда-то украден без осмысления.

что касается указателей на функции то в delphi как-то так это выглядит:

Delphi
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
type
  func = function(const a, b: double): double;
//...
function f(const tau, t: double): double;
begin
  //...
 Result:= /...;
end;
 
function trap(const a, b, eps: double; const afunc: func ): Double;
var
  h1,s,s0,s1,sn, t: double;
  //...
begin
  //...
  s0:=(afunc(a,t)+afunc(b,t))/2;
  //...
  result:= //...;
end;
function DuhamelInt(const double t): double 
const
  eps=0.0000001;  
begin
    // this lines don't needed
    // double a,b,er,eps;
    //double f(double, double),s,trap(double,double,double,double(*) (double,double)) ; 
 
  Result:= trap(0,t,eps,@f); // can be Result:= trap(0, t, eps, f);
end;

USERPC01 · 02.09.2018, 00:37 **[ТС]**

Нужно добавить еще вычислитель производной для начальных условий (См. Попов Основы теории цепей, 1985), чтобы работало по определению (первый фай был посвящен подстановке функции от разных аргументов по шаблону).
Пример интегрирования без параметров f(x)=x*x; a=1, b=5 (y =integr(1;5;f(x))=41.333333), метод Ромберга:

trapint.lpr

Pascal
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
  program trapint;
 uses WinCrt, SysUtils  ;
 
type
 func = function(  x: double ): double;
 
 
  function trap (  a, b, eps: double; const intfunc: func  ): Double;
var
  h1,s,s0,s1,sn   : double;
  i,n: integer;
begin
 
s:=1;
sn:=101;
n:=4;
 
s0:=(intfunc(a )+intfunc(b ))/2;
 
s1:=intfunc ( ((a+b)/2)  );
while (Abs(s-sn)>eps) do
begin
 sn:=s;
 h1:=(b-a)/n;
 
    i:=0 ;
    while  i<(n/2)   do
    begin
 
    s1:=s1+intfunc( a+(2*i+1)*h1 ) ;  ;
    i:=i+1;
    end;
 s:=h1*(s0+s1);
n:=n*2;
 end ;
Result:= s;
end;
 
 
 
 function xsquare( x: double   ): double ; //подынтегральная функция 
begin
 
Result:=x*x ;
end;
 
 procedure GetIntegral;
 const
  eps=1e-8;
 var y: real;
 begin
    y:= trap(1,5, eps,@xsquare);
    writeln('y=',y);
    readln;
 end;
 
 
 
 
begin
 // main;
  GetIntegral;
end.

Добавлено через 7 минут

Pascal
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
program trapint;
 uses WinCrt, SysUtils  ;
 
type
 
  func = function(  x: double ): double;
  funcparam = function(  x: double; param : double): double;
 
 
 
 
 
function   p  ( tau : double)  : double;
 
begin
// /*dU/dt, t=tau;*/
 //double R=1000;
 //double C=0.00000001;
 //double k=0.25;
Result:= -20000*exp(-200*tau);
 
 
end;
 
function h( t :double) : double;
 var
   R,C,k,Tau : double;
begin
 
  R:=1000;
  C:=0.00000001;
  k:=0.25;
  Tau:=R*C;
 
Result:= k*(1-Exp(-t/Tau));
 
 
//Result:=  0.005-0.0025*Exp(-100*t);
end;
 
 
function  f( tau: double ; t: double  ) : double;
begin
 
 
Result:= p (tau)*h(t-tau);
 
end;
 
 
 
 
 
 
 
 
 
 
  function trap (  a, b, eps: double; const intfunc: func  ): Double;
var
  h1,s,s0,s1,sn   : double;
  i,n: integer;
begin
 
s:=1;
sn:=101;
n:=4;
 
s0:=(intfunc(a )+intfunc(b ))/2;
 
s1:=intfunc ( ((a+b)/2)  );
while (Abs(s-sn)>eps) do
begin
 sn:=s;
 h1:=(b-a)/n;
 
    i:=0 ;
    while  i<(n/2)   do
    begin
 
    s1:=s1+intfunc( a+(2*i+1)*h1 ) ;  ;
    i:=i+1;
    end;
 s:=h1*(s0+s1);
n:=n*2;
 end ;
Result:= s;
end;
 
 
 function trapparam (  a, b,t, eps: double; const intfunc: funcparam  ): Double;
var
  h1,s,s0,s1,sn   : double;
  i,n: integer;
begin
 
s:=1;
sn:=101;
n:=4;
 
s0:=(intfunc(a,t )+intfunc(b,t ))/2;
 
s1:=intfunc ( ((a+b)/2),t  );
while (Abs(s-sn)>eps) do
begin
 sn:=s;
 h1:=(b-a)/n;
 
    i:=0 ;
    while  i<(n/2)   do
    begin
 
    s1:=s1+intfunc( ( a+(2*i+1)*h1 ),t) ;
    i:=i+1;
    end;
 s:=h1*(s0+s1);
n:=n*2;
 end ;
Result:= s;
end;
 
 
 
 function xsquare( x: double   ): double ;
begin
 
Result:=x*x ;
end;
 
 procedure GetIntegral;
 const
  eps=1e-8;
 var y: real;
 begin
    y:= trap(1,5, eps,@xsquare);
    writeln('y=',y);
    readln;
 end;
 
 
 
 
 
 
 
 
 
 function DuhamelInt(const  t: double): double;
const
  eps=0.0000001;
begin
 
 
  Result:={ diffsub(t,eps)+ } trapparam(0,t, t,eps,@f); // can be Result:= trap(0, t, eps, f);
end;
 
 
 
procedure  main  ;
 var
      fp : TextFile;
      time,Tend,step,s: double;
begin
 
writeln (' s(t)= trap(0,t,eps,f);   f(tau)=p(tau)*h(t-tau) ');
writeln ('  Input Tend, ms (120) ');
readln( Tend);
Tend:=Tend*0.001;
writeln (' Input step ,ms,(0.01) ');
readln( step);
 
 step:=step*0.001;
 
     assignfile(fp, 'Integral.txt');
      ReWrite(fp);
 
writeln (fp,' s(t)= trap(0,t,eps,f);   f(tau)=p(tau)*h(t-tau) ');
writeln (fp,'  Tbeg=0 , Tend=', Tend);
writeln (fp,'   step=',step,'s ' );
writeln (fp,'   p(tau)=-20000*exp(-200*tau); ');
writeln (fp,'   h(t)=0.005-0.0025*exp(-100*t);  ');
 
  time:=0 ;
 
while    (time<=Tend) do
begin
 
s:=DuhamelInt(time);
writeln(' t= ',time,' s u(t)=  '  ,s) ;
writeln (fp,' t=',time, 's u(t)=  ' ,s);
time:=time+step  ;
end;
  CloseFile( fp);
readln;
 
end;
 
 
 
 
 
 
 
 
 
 
begin
  writeln('trap  test');
  GetIntegral;
  writeln('trapparam test');
  main;
 
end.

Добавлено через 21 минуту
https://ru.wikipedia.org/wiki/... 0%BB%D1%8F

Добавлено через 1 час 42 минуты
trapintparam.lpr

Pascal
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
program trapintparam;
 uses WinCrt{, SysUtils } ;
 
type
 
  func = function(  x: double ): double;
 
 
  funcparamdint  =  function  ( tau: double ; t : double; eps  : double  ; const  Uin  : func  ; const h : func  ) : double;
 
 
 
 
 function  Derivative ( t : double; eps : double ; const   uinp : func     )   : double ;
begin
    Result:=  ( uinp(t+eps) -  uinp(t ) )/eps;
 
end;
 
 
 
function   deriveU  ( tau : double; eps1 : double ; const  Uin   :   func     )  : double;
 
begin
 Result:= Derivative (tau, 1e-8,    Uin     );
end;
 
 
 
function  f_tointegr( tau: double ; t : double; eps  : double  ;  const Uin  : func   ; const  h : func   ) : double;
 const
   eps1=1e-8;
begin
 Result:= deriveU(tau ,eps,  Uin    )*h(t-tau);
end;
 
 
 
 
 
 
 
 
 function trapparam (  a, b,t, eps: double;  const  intfunc: funcparamdint  ;  const uinp : func   ;  const h  : func   ): double;
var
  h1,s,s0,s1,sn ,fun1,fun2,fun3  : double;
  i,n: integer;
begin
 
s:=1;
sn:=101;
n:=4;
   fun1:=intfunc(a,t ,eps  ,  uinp  ,  h     );
   fun2:=intfunc(b,t ,eps  ,  uinp  ,  h    );
s0:=(fun1+fun2)/2;
 
s1:=intfunc ( ((a+b)/2),t ,eps   ,   uinp  ,  h    );
while (Abs(s-sn)>eps) do
begin
 sn:=s;
 h1:=(b-a)/n;
 
    i:=0 ;
    while  i<(n/2)   do
    begin
       fun3:= intfunc( ( a+(2*i+1)*h1 ) , t  ,eps ,  uinp  , h   ) ;
    s1:=s1+fun3;
    i:=i+1;
    end;
 s:=h1*(s0+s1);
n:=n*2;
 end ;
Result:= s;
end;
     { **************************************** }
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 function DuhamelInt(   t: double ; const Uinp: func; const h:func ; eps  : double  ): double;
 
var
   dUinp_Ht :  double;
begin
              //U(t):=U(0)*h(t)+  integr(0;t;derivU(tau)*h(t-tau); dtau)
//U(t):=U(0)*h(t)+  integr(0;t;derivU(tau)*h(t-tau); dtau)
 {
 
   Y1(t) =U1(0)*h(t)+  integr(0;t;derivU(tau)*h(t-tau); dtau)  //0..t1
 
   Y2(t)=U1(0)*h(t)+integr(0;t1;derivU1(tau)*h(t-tau); dtau)+
   +(U2(t1)-U1(t))*(h(t-t1))  +integr(t1; t ; derivU2(tau)*h(t-tau); dtau)+
 
   Y3= U1(0)*h(t)  +  integr(0;t1;derivU1(tau)*h(t-tau); dtau) +
      +(U2(t1)-U1(t1))*(h(t-t1)) +integr(t1; t2 ; derivU2(tau)*h(t-tau); dtau)+
      +(U3(t2)-U2(t2))*(h(t-t2)) +integr(t2; t  ; derivU3(tau)*h(t-tau); dtau)
 
}
        //for dummy  example , if U(t):=U(0)*h(t)+  integr(0;t;derivU(tau)*h(t-tau); dtau)
 
           dUinp_ht:=Uinp (0)*h(t);
 
           Result:=dUinp_Ht+trapparam(0,t, t,eps, @f_tointegr  ,  Uinp    ,  h   );
 
 //
end;
 
 
 
 
    function h( t :double) : double;
  var
    R,C,k,Tau : double;
 begin
 
   R:=1e+3;
   C:=0.1e-6;
   k:=0.25;
   Tau:=R*C;
 
 Result:= k*(1-Exp(-t/Tau));
 
 end;
        function Uinp( t: double ) : double ;
   begin
 
      Result:=1;
   end;
 
 
 
 
 
 
 
procedure  main  ;
 var
      fp : TextFile;
      time,Tend,step,s: double;
begin
 
writeln (' s(t)= trap(0,t,eps,f);   f(tau)=p(tau)*h(t-tau) ');
writeln ('  Input Tend, ms (10) ');
readln( Tend);
Tend:=Tend*0.001;
writeln (' Input step ,ms,(0.1) ');
readln( step);
 
 step:=step*0.001;
 
     assignfile(fp, 'Integral.txt');
      ReWrite(fp);
 
writeln (fp,' s(t)= trap(0,t,eps,f);   f(tau)=p(tau)*h(t-tau) ');
writeln (fp,'  Tbeg=0 , Tend=', Tend);
writeln (fp,'   step=',step,'s ' );
writeln (fp,'   p(tau)=-20000*exp(-200*tau); ');
writeln (fp,'   h(t)=0.005-0.0025*exp(-100*t);  ');
 
  time:=0 ;
 
while    (time<=Tend) do
begin
 
s:=DuhamelInt(time,@Uinp, @h ,1e-6 );
writeln(' t= ',time,' s u(t)=  '  ,s) ;
writeln (fp,' t=',time, 's u(t)=  ' ,s);
time:=time+step  ;
end;
  CloseFile( fp);
readln;
 
end;
 
 
 
 
 
 
 
 
       {
 
 
 
 
 
  function trap (  a, b, eps: double; const intfunc: func  ): Double;
var
  h1,s,s0,s1,sn   : double;
  i,n: integer;
begin
 
s:=1;
sn:=101;
n:=4;
 
s0:=(intfunc(a )+intfunc(b ))/2;
 
s1:=intfunc ( ((a+b)/2)  );
while (Abs(s-sn)>eps) do
begin
 sn:=s;
 h1:=(b-a)/n;
 
    i:=0 ;
    while  i<(n/2)   do
    begin
 
    s1:=s1+intfunc( a+(2*i+1)*h1 ) ;  ;
    i:=i+1;
    end;
 s:=h1*(s0+s1);
n:=n*2;
 end ;
Result:= s;
end;
 
 
 
 
 
 
 
 
 
 function xsquare( x: double   ): double ;
begin
 
Result:=x*x ;
end;
 
 procedure GetIntegral;
 const
  eps=1e-8;
 var y: real;
 begin
    y:= trap(1,5, eps,@xsquare);
    writeln('y=',y);
    readln;
 end;
 
 
       }
 
 
begin
  writeln('trap  test');
 // GetIntegral;
  writeln('trapparam test');
  main;
 
end.

Добавлено через 16 минут
Можно ввести дополнительный параметр eps2 отдельно для

Pascal
1
2
 intfunc(difarg, t ,eps2  ,  uinp  ,  h     ); 
...

, добавив его в аргументы параметрического интегрирования отдельно от

Pascal
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
trapparam (  a, b,t, eps , eps2: double;  const  intfunc: funcparamdint  ;  const uinp : func   ;  const h  : func   ): double;
var
  h1,s,s0,s1,sn ,fun1,fun2,fun3  : double;
  i,n: integer;
begin
 
s:=1;
sn:=101;
n:=4;
   fun1:=intfunc(a,t ,eps2  ,  uinp  ,  h     );
   fun2:=intfunc(b,t ,eps2  ,  uinp  ,  h    );
s0:=(fun1+fun2)/2;
 
s1:=intfunc ( ((a+b)/2),t ,eps2   ,   uinp  ,  h    );
while (Abs(s-sn)>eps) do
begin
 sn:=s;
 h1:=(b-a)/n;
 
    i:=0 ;
    while  i<(n/2)   do
    begin
       fun3:= intfunc( ( a+(2*i+1)*h1 ) , t  ,eps2 ,  uinp  , h   ) ;
    s1:=s1+fun3;
    i:=i+1;
    end;
 s:=h1*(s0+s1);
n:=n*2;
 end ;
Result:= s;
end;
     { **************************************** }
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 function DuhamelInt(   t: double ;  eps ,eps2 : double;  const Uinp: func; const h:func   ): double;
 
var
   dUinp_Ht :  double;
begin
              //U(t):=U(0)*h(t)+  integr(0;t;derivU(tau)*h(t-tau); dtau)
//U(t):=U(0)*h(t)+  integr(0;t;derivU(tau)*h(t-tau); dtau)
 {
 
   Y1(t) =U1(0)*h(t)+  integr(0;t;derivU(tau)*h(t-tau); dtau)  //0..t1
 
   Y2(t)=U1(0)*h(t)+integr(0;t1;derivU1(tau)*h(t-tau); dtau)+
   +(U2(t1)-U1(t))*(h(t-t1))  +integr(t1; t ; derivU2(tau)*h(t-tau); dtau)+
 
   Y3= U1(0)*h(t)  +  integr(0;t1;derivU1(tau)*h(t-tau); dtau) +
      +(U2(t1)-U1(t1))*(h(t-t1)) +integr(t1; t2 ; derivU2(tau)*h(t-tau); dtau)+
      +(U3(t2)-U2(t2))*(h(t-t2)) +integr(t2; t  ; derivU3(tau)*h(t-tau); dtau)
 
}
        //for dummy  example , if U(t):=U(0)*h(t)+  integr(0;t;derivU(tau)*h(t-tau); dtau)
 
           dUinp_ht:=Uinp (0)*h(t);
 
           Result:=dUinp_Ht+trapparam(0,t, t,eps,eps2, @f_tointegr  ,  Uinp    ,  h   );
 
 //
end;
 
 
 
 
    function h( t :double) : double;
  var
    R,C,k,Tau : double;
 begin
 
   R:=1e+3;
   C:=0.1e-6;
   k:=0.25;
   Tau:=R*C;
 
 Result:= k*(1-Exp(-t/Tau));
 
 end;
        function Uinp( t: double ) : double ;
   begin
 
      Result:=1;
   end;
 
 
 
 
 
 
 
procedure  main  ;
 var
      fp : TextFile;
      time,Tend,step,s: double;
begin
 
writeln (' s(t)= trap(0,t,eps,f);   f(tau)=p(tau)*h(t-tau) ');
writeln ('  Input Tend, ms (10) ');
readln( Tend);
Tend:=Tend*0.001;
writeln (' Input step ,ms,(0.1) ');
readln( step);
 
 step:=step*0.001;
 
     assignfile(fp, 'Integral.txt');
      ReWrite(fp);
 
writeln (fp,' s(t)= trap(0,t,eps,f);   f(tau)=p(tau)*h(t-tau) ');
writeln (fp,'  Tbeg=0 , Tend=', Tend);
writeln (fp,'   step=',step,'s ' );
writeln (fp,'   p(tau)=-20000*exp(-200*tau); ');
writeln (fp,'   h(t)=0.005-0.0025*exp(-100*t);  ');
 
  time:=0 ;
 
while    (time<=Tend) do
begin
 
s:=DuhamelInt(time,  1e-6 ,1e-8, @Uinp, @h);
writeln(' t= ',time,' s u(t)=  '  ,s) ;
writeln (fp,' t=',time, 's u(t)=  ' ,s);
time:=time+step  ;
end;
  CloseFile( fp);
readln;
 
end;

Добавлено через 4 минуты
Можно еще ввести массив функций Uinp[i] а также h[j] как объектов , ввести массив вычисленных dUinp_htх[i], чтобы рассчитывать реакции на сложные воздействия , аналогичные

Code
1
2
3
4
5
6
7
8
9
10
11
12
{
 
   Y1(t) =U1(0)*h(t)+  integr(0;t;derivU(tau)*h(t-tau); dtau)  //0..t1
 
   Y2(t)=U1(0)*h(t)+integr(0;t1;derivU1(tau)*h(t-tau); dtau)+
   +(U2(t1)-U1(t))*(h(t-t1))  +integr(t1; t ; derivU2(tau)*h(t-tau); dtau)+
 
   Y3= U1(0)*h(t)  +  integr(0;t1;derivU1(tau)*h(t-tau); dtau) +
      +(U2(t1)-U1(t1))*(h(t-t1)) +integr(t1; t2 ; derivU2(tau)*h(t-tau); dtau)+
      +(U3(t2)-U2(t2))*(h(t-t2)) +integr(t2; t  ; derivU3(tau)*h(t-tau); dtau)
 
}

как в случае https://ru.wikipedia.org/wiki/... 0%BB%D1%8F

Добавлено через 1 час 4 минуты
Для сложных воздействий (при условии типизации объектов, правильной программе инициализации , если устранить ошибки, связанные с типами и поинтерами ) возможна организация следующего рекурсивного механизма вычисления отклика на воздействие для массива функций на интервалах

Pascal
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
 
 function DuhamelInt(   t: double ;  eps , eps2 : double ; const Uinp: func; const h:func  ): double;
type
  Signalarray=record
        Uinp : func;
 
  end;
 
var
   dUinp_Ht :  double;
   y: double;
   tn : array [0..0] of double;
   Signal: Signalarray ;
   n : integer;
 
 
begin
              //U(t):=U(0)*h(t)+  integr(0;t;derivU(tau)*h(t-tau); dtau)
//U(t):=U(0)*h(t)+  integr(0;t;derivU(tau)*h(t-tau); dtau)
 {
 
   Y1(t) =U1(0)*h(t)+  integr(0;t;derivU(tau)*h(t-tau); dtau)  //0..t1
 
   Y2(t)=U1(0)*h(t)+integr(0;t1;derivU1(tau)*h(t-tau); dtau)+
   +(U2(t1)-U1(t))*(h(t-t1))  +integr(t1; t ; derivU2(tau)*h(t-tau); dtau)+
 
   Y3= U1(0)*h(t)  +  integr(0;t1;derivU1(tau)*h(t-tau); dtau) +
      +(U2(t1)-U1(t1))*(h(t-t1)) +integr(t1; t2 ; derivU2(tau)*h(t-tau); dtau)+
      +(U3(t2)-U2(t2))*(h(t-t2)) +integr(t2; t  ; derivU3(tau)*h(t-tau); dtau)
 
}
        //for dummy  example , if U(t):=U(0)*h(t)+  integr(0;t;derivU(tau)*h(t-tau); dtau)
 
           n:=3;
 
         SetLength(tn,n+1) ;    // 0,1,2,3  t[0]=0 t[1]=t1,t[2]=t2, t[3]=t
         SetLength(Signal,n);   //0,1,2- >U1,U2,U3
 
         //  ReadSignalArray(tn,Signal); ?
 
 
             {
 
 
 
           tn[0]:=0 ;
           Signal.Uinp [0]:=U1;       //object, function U1(t) правильно присвоить тип, ассоциировать 
 
 
 
           tn[1]:=t1 ;      //t1       точки разрыва
           Signal.Uinp [1]:=U2;          //object, function U2(t) 
 
 
           tn[2]:=t2 ;        //t2
           Signal.Uinp [2]:=U3;             //    //object, function U3(t), ...
 
 
            tn[n]:=t ;
 
            }
 
 
            y:=0;
            dUinp_ht := Signal[0].Uinp(tn[0]) *CircuitParam.h(t-tn[0] );
            y := dUinp_ht +trapparam(tn[0],tn[1], t,eps,eps2, @f_tointegr  ,  Uinput[0]    ,  h   );
 
 
 
           for i:=1 to n-1 do
            begin
 
              dUinp_ht[i]:=( Signal[i].Uinp(tn[i])-Signal[i].Uinp (tn[i-1])  )*CircuitParam.h(t-tn[i]);
              y :=y +dUinp_ht[i]+trapparam(tn[i],tn[i+1], t,eps,eps2, @f_tointegr  , Uinput[i]    ,  h   );
 
            end;
 
           Result:=y;
 
 
          // Result:=dUinp_Ht+trapparam(0,t, t,eps,eps2, @f_tointegr  ,  Uinp    ,  h   );
 
 
end;

Добавлено через 1 минуту
Как организовать этот алгоритм подстановки функций корректно используя в том числе раздел типов, корректную передачу адресов подпрограмм вычисления вектора воздействия (как объектов)?

Добавлено через 2 часа 24 минуты
Вариант с несколькими сигналами (возможно , с багами (t>tmax4; корректно ли считает в точках разрыва ), может, что-то в принципе действия и в точках разрыва исправить , на уровне формул перепроверить ?). Объекты размещаемы в памяти .

Pascal
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
program paramtrapint;
 uses WinCrt{, SysUtils}  ;
 
type
 
  func = function(  x: double ): double;
   Signalarray=record
        Uinp : func;
             tbegin : double;
             tend: double;
     end;
 
 
  funcparamdint  =  function  ( tau: double ; t : double; eps  : double  ; const  Uin  : func  ; const h : func  ) : double;
 
 
 
 
 function  Derivative ( t : double; eps : double ; const   uinp : func     )   : double ;
begin
    Result:=  ( uinp(t+eps) -  uinp(t ) )/eps;
 
end;
 
 
 
function   deriveU  ( tau : double; eps1 : double ; const  Uin   :   func     )  : double;
 
begin
 Result:= Derivative (tau, 1e-8,    Uin     );
end;
 
 
 
function  f_tointegr( tau: double ; t : double; eps  : double  ;  const Uin  : func   ; const  h : func   ) : double;
 const
   eps1=1e-8;
begin
 Result:= deriveU(tau ,eps,  Uin    )*h(t-tau);
end;
 
 
 
 
 
 
 
 
 function trapparam (  a, b,t, eps, eps2: double;  const  intfunc: funcparamdint  ;  const uinp : func   ;  const h  : func   ): double;
var
  h1,s,s0,s1,sn ,fun1,fun2,fun3  : double;
  i,n: integer;
begin
 
s:=1;
sn:=101;
n:=4;
   fun1:=intfunc(a,t ,eps2  ,  uinp  ,  h     );
   fun2:=intfunc(b,t ,eps2  ,  uinp  ,  h    );
s0:=(fun1+fun2)/2;
 
s1:=intfunc ( ((a+b)/2),t ,eps2   ,   uinp  ,  h    );
while (Abs(s-sn)>eps) do
begin
 sn:=s;
 h1:=(b-a)/n;
 
    i:=0 ;
    while  i<(n/2)   do
    begin
       fun3:= intfunc( ( a+(2*i+1)*h1 ) , t  ,eps2 ,  uinp  , h   ) ;
    s1:=s1+fun3;
    i:=i+1;
    end;
 s:=h1*(s0+s1);
n:=n*2;
 end ;
Result:= s;
end;
     { **************************************** }
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
      function h( t :double) : double;
  var
    R,C,k,Tau : double;
 begin
 
   R:=1e+3;
   C:=0.1e-6;
   k:=1;   //0.25
   Tau:=R*C;
   if(t>=0)  then Result:= (k*(1-Exp(-t/Tau)))  else Result:=0;
 
 end;
 
 
 
    function Uinp1( t: double ) : double ;
   begin
 
      Result:=1 ;  //t 0...0.1 s
   end;
 
     function Uinp2( t: double ) : double ;
   begin
 
      Result:=0;   // 0.1 s...3 s
   end;
 
 
       function Uinp3( t: double ) : double ;
   begin
 
      Result:=-1; // 3s..3.2 s
   end;
 
     function Uinp4( t: double ) : double ;
   begin
 
      Result:=0;   // 3.2s..6 s
   end;
 
   function AssignSignal(const   Ui  : func; tbeg :double; tend: double ): Signalarray;
   var
     Uinput :Signalarray ;
   begin
 
              Uinput.Uinp:=Ui;
              Uinput.tbegin:=tbeg;
              Uinput.tend:=tend;
 
      Result:=Uinput;
   end;
 
 
 
  function getht( const h: func ; t : double)  : double;
begin
       Result:=h(t);
end;
 
 
 
 
 function DuhamelInt(   t: double ;  eps , eps2 : double  ;   const h:func ): double;
 
 
 
var
   dUinp_Ht :  double;
   y: double;
   te  ,tp : double;
   Signal: array of Signalarray ;
   n,m ,i : integer;
 
 
begin
              //U(t):=U(0)*h(t)+  integr(0;t;derivU(tau)*h(t-tau); dtau)
//U(t):=U(0)*h(t)+  integr(0;t;derivU(tau)*h(t-tau); dtau)
 {
    u
 
 
    t=[0 ,t1]
   Y1(t) =U1(0)*h(t)+  integr(0;t;derivU(tau)*h(t-tau); dtau)  //0..t1
     Отклик на остальных интервалах вычисляется по формулам, вытекающих из принципа суперпозиции:
   Y2(t)=U1(0)*h(t)+integr(0;t1;derivU1(tau)*h(t-tau); dtau)+
   +(U2(t1)-U1(t1))*(h(t-t1))  +integr(t1; t ; derivU2(tau)*h(t-tau); dtau)+
 
   Y3= U1(0)*h(t)  +  integr(0;t1;derivU1(tau)*h(t-tau); dtau) +
      +(U2(t1)-U1(t1))*(h(t-t1)) +integr(t1; t2 ; derivU2(tau)*h(t-tau); dtau)+
      +(U3(t2)-U2(t2))*(h(t-t2)) +integr(t2; t  ; derivU3(tau)*h(t-tau); dtau)
 
}
        //for dummy  example , if U(t):=U(0)*h(t)+  integr(0;t;derivU(tau)*h(t-tau); dtau)
 
           n:=4;
           //n:=4
 
      //   SetLength(tn,n+1) ;    // 0,1,2,3  t[0]=0 t[1]=t1,t[2]=t2, t[3]=t
         SetLength(Signal,n);   //0,1,2- >U1,U2,U3
 
           Signal[0]:=AssignSignal(  @Uinp1  ,0, 1e-3  ) ;   //0,t1
           Signal[1]:=AssignSignal(  @Uinp2  ,1e-3,2e-3   ) ;   //t1,t2
           Signal[2]:=AssignSignal(  @Uinp3  ,2e-3,3   ) ;   //t2,t3
           Signal[3]:=AssignSignal(  @Uinp4  ,3e-3,4e-3   ) ;   //t3,t4
 
 
 
 
 
          m:=0;
         if(t>=0   ) and   (t<= Signal[0].tend) then m:=0;
         if(t>= Signal[1].tbegin) and   (t<= Signal[1].tend) then m:=1;
         if(t>= Signal[2].tbegin) and   (t<= Signal[2].tend) then m:=2;
         if(t>= Signal[3].tbegin) and   (t<= Signal[3].tend) then m:=3;
 
            y:=0;
 
           for i:=0 to m do
            begin
 
 
 
 
 
 
                  if (i=0) then
                 begin
                     tp:=0  ;
                     te:=t  ;
                     dUinp_ht := Signal[0].Uinp(tp) * getht(h,t-tp)    ;
                    y := dUinp_ht +trapparam(tp,te,t,eps,eps2, @f_tointegr,Signal[0].Uinp,  h    );
                 end
 
                  else
                   begin
                          tp:=Signal[i].tbegin  ;  //t1,t2 ;
                          te:=Signal[i].tend  ;    //t2,t3 ;
                          if (i=m)  then   te:=t;
                          dUinp_ht:=( Signal[i].Uinp(tp)-Signal[i-1].Uinp (tp))*getht(h,t-tp);
                          y :=y+dUinp_ht+trapparam(tp,te, t,eps,eps2, @f_tointegr,Signal[i].Uinp , h  );
                     end;
            end;
 
           Result:=y;
 
 
          // Result:=dUinp_Ht+trapparam(0,t, t,eps,eps2, @f_tointegr  ,  Uinp    ,  h   );
 
 
end;
 
 
 
 
 
 
 
 
 
 
 
 
procedure  main  ;
 var
      fp : TextFile;
      time,Tend,step,s: double;
begin
 
writeln (' s(t)= trap(0,t,eps,f);   f(tau)=p(tau)*h(t-tau) ');
writeln ('  Input Tend, ms (10) ');
readln( Tend);
Tend:=Tend*0.001 ;
writeln (' Input step ,ms,(0.1) ');
readln( step);
 
 step:=step*0.001;
 
     assignfile(fp, 'Integral.txt');
      ReWrite(fp);
 
writeln (fp,' s(t)= trap(0,t,eps,f);   f(tau)=p(tau)*h(t-tau) ');
writeln (fp,'  Tbeg=0 , Tend=', Tend);
writeln (fp,'   step=',step,'s ' );
writeln (fp,'   p(tau)=-20000*exp(-200*tau); ');
writeln (fp,'   h(t)=0.005-0.0025*exp(-100*t);  ');
 
  time:=0 ;
 
while    (time<=Tend) do
begin
 
s:=DuhamelInt(time, 1e-6,1e-8,   @h   );
writeln(' t= ',time,' s u(t)=  '  ,s) ;
writeln (fp,' t=',time, 's u(t)=  ' ,s);
time:=time+step  ;
end;
  CloseFile( fp);
readln;
 
end;
 
 
 
 
 
 
 
 
       {
 
 
 
 
 
  function trap (  a, b, eps: double; const intfunc: func  ): Double;
var
  h1,s,s0,s1,sn   : double;
  i,n: integer;
begin
 
s:=1;
sn:=101;
n:=4;
 
s0:=(intfunc(a )+intfunc(b ))/2;
 
s1:=intfunc ( ((a+b)/2)  );
while (Abs(s-sn)>eps) do
begin
 sn:=s;
 h1:=(b-a)/n;
 
    i:=0 ;
    while  i<(n/2)   do
    begin
 
    s1:=s1+intfunc( a+(2*i+1)*h1 ) ;  ;
    i:=i+1;
    end;
 s:=h1*(s0+s1);
n:=n*2;
 end ;
Result:= s;
end;
 
 
 
 
 
 
 
 
 
 function xsquare( x: double   ): double ;
begin
 
Result:=x*x ;
end;
 
 procedure GetIntegral;
 const
  eps=1e-8;
 var y: real;
 begin
    y:= trap(1,5, eps,@xsquare);
    writeln('y=',y);
    readln;
 end;
 
 
       }
 
 
begin
  writeln('trap  test');
 // GetIntegral;
  writeln('trapparam test');
  main;
 
  readln;
 
end.

USERPC01 · 02.09.2018, 01:21 **[ТС]**

результаты

Code
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
trap  test
trapparam test
 s(t)= trap(0,t,eps,f);   f(tau)=p(tau)*h(t-tau)
  Input Tend, ms (10)
10
 Input step ,ms,(0.1)
0.1
 t=  0.0000000000000000E+000 s u(t)=   0.0000000000000000E+000
 t=  1.0000000000000000E-004 s u(t)=   6.3212055882855778E-001
 t=  2.0000000000000001E-004 s u(t)=   8.6466471676338730E-001
 t=  3.0000000000000003E-004 s u(t)=   9.5021293163213605E-001
 t=  4.0000000000000002E-004 s u(t)=   9.8168436111126578E-001
 t=  5.0000000000000001E-004 s u(t)=   9.9326205300091452E-001
 t=  6.0000000000000006E-004 s u(t)=   9.9752124782333362E-001
 t=  7.0000000000000010E-004 s u(t)=   9.9908811803444553E-001
 t=  8.0000000000000015E-004 s u(t)=   9.9966453737209748E-001
 t=  9.0000000000000019E-004 s u(t)=   9.9987659019591335E-001
 t=  1.0000000000000002E-003 s u(t)=   9.9995460007023529E-001
 t=  1.1000000000000003E-003 s u(t)=   3.6786273947065107E-001
 t=  1.2000000000000003E-003 s u(t)=   1.3532913902425892E-001
 t=  1.3000000000000004E-003 s u(t)=   4.9784808038456752E-002
 t=  1.4000000000000004E-003 s u(t)=   1.8314807360015006E-002
 t=  1.5000000000000005E-003 s u(t)=   6.7376410967650013E-003
 t=  1.6000000000000005E-003 s u(t)=   2.4786396414916423E-003
 t=  1.7000000000000006E-003 s u(t)=   9.1184056617732434E-004
 t=  1.8000000000000006E-003 s u(t)=   3.3544739792279454E-004
 t=  1.9000000000000006E-003 s u(t)=   1.2340420129020035E-004
 t=  2.0000000000000005E-003 s u(t)=   4.5397868604572789E-005
 t=  2.1000000000000003E-003 s u(t)=  -6.3210385788602441E-001
 t=  2.2000000000000001E-003 s u(t)=  -8.6465857282998082E-001
 t=  2.3000000000000000E-003 s u(t)=  -9.5021067140534787E-001
 t=  2.3999999999999998E-003 s u(t)=  -9.8168352962029803E-001
 t=  2.4999999999999996E-003 s u(t)=  -9.9326174711248194E-001
 t=  2.5999999999999994E-003 s u(t)=  -9.9752113529326802E-001
 t=  2.6999999999999993E-003 s u(t)=  -9.9908807663694787E-001
 t=  2.7999999999999991E-003 s u(t)=  -9.9966452214280921E-001
 t=  2.8999999999999989E-003 s u(t)=  -9.9987658459337125E-001
 t=  2.9999999999999988E-003 s u(t)=  -9.9995459800917752E-001
 t=  3.0999999999999986E-003 s u(t)=  -3.6786273871243580E-001
 t=  3.1999999999999984E-003 s u(t)=  -1.3532913874532748E-001
 t=  3.2999999999999982E-003 s u(t)=  -4.9784807935843722E-002
 t=  3.3999999999999981E-003 s u(t)=  -1.8314807322265758E-002
 t=  3.4999999999999979E-003 s u(t)=  -6.7376410828778877E-003
 t=  3.5999999999999977E-003 s u(t)=  -2.4786396363827290E-003
 t=  3.6999999999999976E-003 s u(t)=  -9.1184056429804983E-004
 t=  3.7999999999999974E-003 s u(t)=  -3.3544739723134764E-004
 t=  3.8999999999999972E-003 s u(t)=  -1.2340420103584826E-004
 t=  3.9999999999999975E-003 s u(t)=  -4.5397868515317796E-005
 t=  4.0999999999999977E-003 s u(t)=  -9.9999999924170957E-001
 t=  4.1999999999999980E-003 s u(t)=  -9.9999999972104059E-001
 t=  4.2999999999999983E-003 s u(t)=  -9.9999999989737653E-001
 t=  4.3999999999999985E-003 s u(t)=  -9.9999999996224698E-001
 t=  4.4999999999999988E-003 s u(t)=  -9.9999999998611144E-001
 t=  4.5999999999999991E-003 s u(t)=  -9.9999999999489064E-001
 t=  4.6999999999999993E-003 s u(t)=  -9.9999999999812039E-001
 t=  4.7999999999999996E-003 s u(t)=  -9.9999999999930855E-001
 t=  4.8999999999999998E-003 s u(t)=  -9.9999999999974565E-001
 t=  5.0000000000000001E-003 s u(t)=  -9.9999999999990641E-001
 t=  5.1000000000000004E-003 s u(t)=  -9.9999999999996558E-001
 t=  5.2000000000000006E-003 s u(t)=  -9.9999999999998734E-001
 t=  5.3000000000000009E-003 s u(t)=  -9.9999999999999534E-001
 t=  5.4000000000000012E-003 s u(t)=  -9.9999999999999833E-001
 t=  5.5000000000000014E-003 s u(t)=  -9.9999999999999933E-001
 t=  5.6000000000000017E-003 s u(t)=  -9.9999999999999978E-001
 t=  5.7000000000000019E-003 s u(t)=  -9.9999999999999989E-001
 t=  5.8000000000000022E-003 s u(t)=  -1.0000000000000000E+000
 t=  5.9000000000000025E-003 s u(t)=  -1.0000000000000000E+000
 t=  6.0000000000000027E-003 s u(t)=  -1.0000000000000000E+000
 t=  6.1000000000000030E-003 s u(t)=  -1.0000000000000000E+000
 t=  6.2000000000000033E-003 s u(t)=  -1.0000000000000000E+000
 t=  6.3000000000000035E-003 s u(t)=  -1.0000000000000000E+000
 t=  6.4000000000000038E-003 s u(t)=  -1.0000000000000000E+000
 t=  6.5000000000000040E-003 s u(t)=  -1.0000000000000000E+000
 t=  6.6000000000000043E-003 s u(t)=  -1.0000000000000000E+000
 t=  6.7000000000000046E-003 s u(t)=  -1.0000000000000000E+000
 t=  6.8000000000000048E-003 s u(t)=  -1.0000000000000000E+000
 t=  6.9000000000000051E-003 s u(t)=  -1.0000000000000000E+000
 t=  7.0000000000000053E-003 s u(t)=  -1.0000000000000000E+000
 t=  7.1000000000000056E-003 s u(t)=  -1.0000000000000000E+000
 t=  7.2000000000000059E-003 s u(t)=  -1.0000000000000000E+000
 t=  7.3000000000000061E-003 s u(t)=  -1.0000000000000000E+000
 t=  7.4000000000000064E-003 s u(t)=  -1.0000000000000000E+000
 t=  7.5000000000000067E-003 s u(t)=  -1.0000000000000000E+000
 t=  7.6000000000000069E-003 s u(t)=  -1.0000000000000000E+000
 t=  7.7000000000000072E-003 s u(t)=  -1.0000000000000000E+000
 t=  7.8000000000000074E-003 s u(t)=  -1.0000000000000000E+000
 t=  7.9000000000000077E-003 s u(t)=  -1.0000000000000000E+000
 t=  8.0000000000000071E-003 s u(t)=  -1.0000000000000000E+000
 t=  8.1000000000000065E-003 s u(t)=  -1.0000000000000000E+000
 t=  8.2000000000000059E-003 s u(t)=  -1.0000000000000000E+000
 t=  8.3000000000000053E-003 s u(t)=  -1.0000000000000000E+000
 t=  8.4000000000000047E-003 s u(t)=  -1.0000000000000000E+000
 t=  8.5000000000000041E-003 s u(t)=  -1.0000000000000000E+000
 t=  8.6000000000000035E-003 s u(t)=  -1.0000000000000000E+000
 t=  8.7000000000000029E-003 s u(t)=  -1.0000000000000000E+000
 t=  8.8000000000000023E-003 s u(t)=  -1.0000000000000000E+000
 t=  8.9000000000000017E-003 s u(t)=  -1.0000000000000000E+000
 t=  9.0000000000000011E-003 s u(t)=  -1.0000000000000000E+000
 t=  9.1000000000000004E-003 s u(t)=  -1.0000000000000000E+000
 t=  9.1999999999999998E-003 s u(t)=  -1.0000000000000000E+000
 t=  9.2999999999999992E-003 s u(t)=  -1.0000000000000000E+000
 t=  9.3999999999999986E-003 s u(t)=  -1.0000000000000000E+000
 t=  9.4999999999999980E-003 s u(t)=  -1.0000000000000000E+000
 t=  9.5999999999999974E-003 s u(t)=  -1.0000000000000000E+000
 t=  9.6999999999999968E-003 s u(t)=  -1.0000000000000000E+000
 t=  9.7999999999999962E-003 s u(t)=  -1.0000000000000000E+000
 t=  9.8999999999999956E-003 s u(t)=  -1.0000000000000000E+000
 t=  9.9999999999999950E-003 s u(t)=  -1.0000000000000000E+000
^A

Добавлено через 32 минуты
Как уменьшить ошибку шага с t= 3.0000000000000003E-004 s (если не использовать форматированный вывод , double мало что ли, с type double=real ,uses Math, uses SysUtils, uses Crt ,uses WinCrt то же ) ? Урезать можно, но хранятся ли данные без ошибок ? Обычно double хватало, это с float или real до 5...6 знаков после запятой , особенно в Си++.

Pascal
1
2
3
4
5
6
7
8
 t=  0.0000000000000000E+000 s u(t)=   0.0000000000000000E+000
 t=  1.0000000000000000E-004 s u(t)=   6.3212055882855778E-001
 t=  2.0000000000000001E-004 s u(t)=   8.6466471676338730E-001
 t=  3.0000000000000003E-004 s u(t)=   9.5021293163213605E-001
 t=  4.0000000000000002E-004 s u(t)=   9.8168436111126578E-001
 t=  5.0000000000000001E-004 s u(t)=   9.9326205300091452E-001
 t=  6.0000000000000006E-004 s u(t)=   9.9752124782333362E-001
 t=  7.0000000000000010E-004 s u(t)=   9.9908811803444553E-001

Добавлено через 22 секунды
А в дельфи лучше?

Добавлено через 43 секунды

Code
1
2
3
 trap  test
trapparam test
 s(t)= trap(0,t,eps,f);   f(tau)=p(tau)*h(t-tau)

можно убрать

Добавлено через 2 минуты
С другими постоянными времени

Pascal
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
program paramtrapint;
 uses    SysUtils  ;
 
type
  //double=real;
  func = function(  x: double ): double;
   Signalarray=record
        Uinp : func;
             tbegin : double;
             tend: double;
     end;
 
 
  funcparamdint  =  function  ( tau: double ; t : double; eps  : double  ; const  Uin  : func  ; const h : func  ) : double;
 
 
 
 
 function  Derivative ( t : double; eps : double ; const   uinp : func     )   : double ;
begin
    Result:=  ( uinp(t+eps) -  uinp(t ) )/eps;
 
end;
 
 
 
function   deriveU  ( tau : double; eps1 : double ; const  Uin   :   func     )  : double;
 
begin
 Result:= Derivative (tau, 1e-8,    Uin     );
end;
 
 
 
function  f_tointegr( tau: double ; t : double; eps  : double  ;  const Uin  : func   ; const  h : func   ) : double;
 const
   eps1=1e-8;
begin
 Result:= deriveU(tau ,eps,  Uin    )*h(t-tau);
end;
 
 
 
 
 
 
 
 
 function trapparam (  a, b,t, eps, eps2: double;  const  intfunc: funcparamdint  ;  const uinp : func   ;  const h  : func   ): double;
var
  h1,s,s0,s1,sn ,fun1,fun2,fun3  : double;
  i,n: integer;
begin
 
s:=1;
sn:=101;
n:=4;
   fun1:=intfunc(a,t ,eps2  ,  uinp  ,  h     );
   fun2:=intfunc(b,t ,eps2  ,  uinp  ,  h    );
s0:=(fun1+fun2)/2;
 
s1:=intfunc ( ((a+b)/2),t ,eps2   ,   uinp  ,  h    );
while (Abs(s-sn)>eps) do
begin
 sn:=s;
 h1:=(b-a)/n;
 
    i:=0 ;
    while  i<(n/2)   do
    begin
       fun3:= intfunc( ( a+(2*i+1)*h1 ) , t  ,eps2 ,  uinp  , h   ) ;
    s1:=s1+fun3;
    i:=i+1;
    end;
 s:=h1*(s0+s1);
n:=n*2;
 end ;
Result:= s;
end;
     { **************************************** }
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
      function h( t :double) : double;
  var
    R,C,k,Tau : double;
 begin
 
   R:=100;
   C:=0.1e-6;
   k:=1;   //0.25
   Tau:=R*C;
   if(t>=0)  then Result:= (k*(1-Exp(-t/Tau)))  else Result:=0;
 
 end;
 
 
 
    function Uinp1( t: double ) : double ;
   begin
 
      Result:=1 ;  //t 0...0.1 s
   end;
 
     function Uinp2( t: double ) : double ;
   begin
 
      Result:=0;   // 0.1 s...3 s
   end;
 
 
       function Uinp3( t: double ) : double ;
   begin
 
      Result:=-1; // 3s..3.2 s
   end;
 
     function Uinp4( t: double ) : double ;
   begin
 
      Result:=0;   // 3.2s..6 s
   end;
 
   function AssignSignal(const   Ui  : func; tbeg :double; tend: double ): Signalarray;
   var
     Uinput :Signalarray ;
   begin
 
              Uinput.Uinp:=Ui;
              Uinput.tbegin:=tbeg;
              Uinput.tend:=tend;
 
      Result:=Uinput;
   end;
 
 
 
  function getht( const h: func ; t : double)  : double;
begin
       Result:=h(t);
end;
 
 
 
 
 function DuhamelInt(   t: double ;  eps , eps2 : double  ;   const h:func ): double;
 
 
 
var
   dUinp_Ht :  double;
   y: double;
   te  ,tp : double;
   Signal: array of Signalarray ;
   n,m ,i : integer;
 
 
begin
              //U(t):=U(0)*h(t)+  integr(0;t;derivU(tau)*h(t-tau); dtau)
//U(t):=U(0)*h(t)+  integr(0;t;derivU(tau)*h(t-tau); dtau)
 {
    u
 
 
    t=[0 ,t1]
   Y1(t) =U1(0)*h(t)+  integr(0;t;derivU(tau)*h(t-tau); dtau)  //0..t1
     Отклик на остальных интервалах вычисляется по формулам, вытекающих из принципа суперпозиции:
   Y2(t)=U1(0)*h(t)+integr(0;t1;derivU1(tau)*h(t-tau); dtau)+
   +(U2(t1)-U1(t1))*(h(t-t1))  +integr(t1; t ; derivU2(tau)*h(t-tau); dtau)+
 
   Y3= U1(0)*h(t)  +  integr(0;t1;derivU1(tau)*h(t-tau); dtau) +
      +(U2(t1)-U1(t1))*(h(t-t1)) +integr(t1; t2 ; derivU2(tau)*h(t-tau); dtau)+
      +(U3(t2)-U2(t2))*(h(t-t2)) +integr(t2; t  ; derivU3(tau)*h(t-tau); dtau)
 
}
        //for dummy  example , if U(t):=U(0)*h(t)+  integr(0;t;derivU(tau)*h(t-tau); dtau)
 
           n:=4;
           //n:=4
 
      //   SetLength(tn,n+1) ;    // 0,1,2,3  t[0]=0 t[1]=t1,t[2]=t2, t[3]=t
         SetLength(Signal,n);   //0,1,2- >U1,U2,U3
 
           Signal[0]:=AssignSignal(  @Uinp1  ,0, 1e-3  ) ;   //0,t1
           Signal[1]:=AssignSignal(  @Uinp2  ,1e-3,2e-3   ) ;   //t1,t2
           Signal[2]:=AssignSignal(  @Uinp3  ,2e-3,3   ) ;   //t2,t3
           Signal[3]:=AssignSignal(  @Uinp4  ,3e-3,4e-3   ) ;   //t3,t4
 
 
 
 
 
          m:=0;
         if(t>=0   ) and   (t<= Signal[0].tend) then m:=0;
         if(t>= Signal[1].tbegin) and   (t<= Signal[1].tend) then m:=1;
         if(t>= Signal[2].tbegin) and   (t<= Signal[2].tend) then m:=2;
         if(t>= Signal[3].tbegin) and   (t<= Signal[3].tend) then m:=3;
 
            y:=0;
 
           for i:=0 to m do
            begin
 
 
 
 
 
 
                  if (i=0) then
                 begin
                     tp:=0  ;
                     te:=t  ;
                     dUinp_ht := Signal[0].Uinp(tp) * getht(h,t-tp)    ;
                    y := dUinp_ht +trapparam(tp,te,t,eps,eps2, @f_tointegr,Signal[0].Uinp,  h    );
                 end
 
                  else
                   begin
                          tp:=Signal[i].tbegin  ;  //t1,t2 ;
                          te:=Signal[i].tend  ;    //t2,t3 ;
                          if (i=m)  then   te:=t;
                          dUinp_ht:=( Signal[i].Uinp(tp)-Signal[i-1].Uinp (tp))*getht(h,t-tp);
                          y :=y+dUinp_ht+trapparam(tp,te, t,eps,eps2, @f_tointegr,Signal[i].Uinp , h  );
                     end;
            end;
 
           Result:=y;
 
 
          // Result:=dUinp_Ht+trapparam(0,t, t,eps,eps2, @f_tointegr  ,  Uinp    ,  h   );
 
 
end;
 
 
 
 
 
 
 
 
 
 
 
 
procedure  main  ;
 var
      fp : TextFile;
      time,Tend,step,s: double;
begin
 
//writeln (' s(t)= trap(0,t,eps,f);   f(tau)=p(tau)*h(t-tau) ');
writeln ('  Input Tend, s (0.01) ');
readln( Tend);
//Tend:=Tend*0.001 ;
writeln (' Input step ,s,(0.0001) ');
readln( step);
 
 //step:=step*0.001;
 
     assignfile(fp, 'Integral.txt');
      ReWrite(fp);
 
 
writeln (fp,'  Tbeg=0 , Tend=', Tend);
writeln (fp,'   step=',step,'s ' );
//writeln (fp,'   p(tau)=-20000*exp(-200*tau); ');
//writeln (fp,'   h(t)=0.005-0.0025*exp(-100*t);  ');
 
  time:=0 ;
 
while    (time<=Tend) do
begin
 
s:=DuhamelInt(time, 1e-6,1e-8,   @h   );
writeln(' t= ',time,' s u(t)=  '  ,s) ;
writeln (fp,' t=',time, 's u(t)=  ' ,s);
time:=time+step  ;
end;
  CloseFile( fp);
readln;
 
end;
 
 
 
 
 
 
 
 
       {
 
 
 
 
 
  function trap (  a, b, eps: double; const intfunc: func  ): Double;
var
  h1,s,s0,s1,sn   : double;
  i,n: integer;
begin
 
s:=1;
sn:=101;
n:=4;
 
s0:=(intfunc(a )+intfunc(b ))/2;
 
s1:=intfunc ( ((a+b)/2)  );
while (Abs(s-sn)>eps) do
begin
 sn:=s;
 h1:=(b-a)/n;
 
    i:=0 ;
    while  i<(n/2)   do
    begin
 
    s1:=s1+intfunc( a+(2*i+1)*h1 ) ;  ;
    i:=i+1;
    end;
 s:=h1*(s0+s1);
n:=n*2;
 end ;
Result:= s;
end;
 
 
 
 
 
 
 
 
 
 function xsquare( x: double   ): double ;
begin
 
Result:=x*x ;
end;
 
 procedure GetIntegral;
 const
  eps=1e-8;
 var y: real;
 begin
    y:= trap(1,5, eps,@xsquare);
    writeln('y=',y);
    readln;
 end;
 
 
       }
 
 
begin
  writeln('trap  test');
 // GetIntegral;
  writeln('trapparam test');
  main;
 
  readln;
 
end.

Code
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
trap  test
trapparam test
  Input Tend, s (0.01)
0.01
 Input step ,s,(0.0001)
0.0001
 t=  0.0000000000000000E+000 s u(t)=   0.0000000000000000E+000
 t=  1.0000000000000000E-004 s u(t)=   9.9995460007023751E-001
 t=  2.0000000000000001E-004 s u(t)=   9.9999999793884642E-001
 t=  3.0000000000000003E-004 s u(t)=   9.9999999999990641E-001
 t=  4.0000000000000002E-004 s u(t)=   1.0000000000000000E+000
 t=  5.0000000000000001E-004 s u(t)=   1.0000000000000000E+000
 t=  6.0000000000000006E-004 s u(t)=   1.0000000000000000E+000
 t=  7.0000000000000010E-004 s u(t)=   1.0000000000000000E+000
 t=  8.0000000000000015E-004 s u(t)=   1.0000000000000000E+000
 t=  9.0000000000000019E-004 s u(t)=   1.0000000000000000E+000
 t=  1.0000000000000002E-003 s u(t)=   9.9999999999997835E-001
 t=  1.1000000000000003E-003 s u(t)=   4.5399929762490743E-005
 t=  1.2000000000000003E-003 s u(t)=   2.0611535811454473E-009
 t=  1.3000000000000004E-003 s u(t)=   9.3591800975900696E-014
 t=  1.4000000000000004E-003 s u(t)=   0.0000000000000000E+000
 t=  1.5000000000000005E-003 s u(t)=   0.0000000000000000E+000
 t=  1.6000000000000005E-003 s u(t)=   0.0000000000000000E+000
 t=  1.7000000000000006E-003 s u(t)=   0.0000000000000000E+000
 t=  1.8000000000000006E-003 s u(t)=   0.0000000000000000E+000
 t=  1.9000000000000006E-003 s u(t)=   0.0000000000000000E+000
 t=  2.0000000000000005E-003 s u(t)=  -4.3368086899420177E-014
 t=  2.1000000000000003E-003 s u(t)=  -9.9995460007023751E-001
 t=  2.2000000000000001E-003 s u(t)=  -9.9999999793884642E-001
 t=  2.3000000000000000E-003 s u(t)=  -9.9999999999990641E-001
 t=  2.3999999999999998E-003 s u(t)=  -1.0000000000000000E+000
 t=  2.4999999999999996E-003 s u(t)=  -1.0000000000000000E+000
 t=  2.5999999999999994E-003 s u(t)=  -1.0000000000000000E+000
 t=  2.6999999999999993E-003 s u(t)=  -1.0000000000000000E+000
 t=  2.7999999999999991E-003 s u(t)=  -1.0000000000000000E+000
 t=  2.8999999999999989E-003 s u(t)=  -1.0000000000000000E+000
 t=  2.9999999999999988E-003 s u(t)=  -1.0000000000000000E+000
 t=  3.0999999999999986E-003 s u(t)=  -4.5399929762490743E-005
 t=  3.1999999999999984E-003 s u(t)=  -2.0611535811454473E-009
 t=  3.2999999999999982E-003 s u(t)=  -9.3591800975900696E-014
 t=  3.3999999999999981E-003 s u(t)=   0.0000000000000000E+000
 t=  3.4999999999999979E-003 s u(t)=   0.0000000000000000E+000
 t=  3.5999999999999977E-003 s u(t)=   0.0000000000000000E+000
 t=  3.6999999999999976E-003 s u(t)=   0.0000000000000000E+000
 t=  3.7999999999999974E-003 s u(t)=   0.0000000000000000E+000
 t=  3.8999999999999972E-003 s u(t)=   0.0000000000000000E+000
 t=  3.9999999999999975E-003 s u(t)=   0.0000000000000000E+000
 t=  4.0999999999999977E-003 s u(t)=  -1.0000000000000000E+000
 t=  4.1999999999999980E-003 s u(t)=  -1.0000000000000000E+000
 t=  4.2999999999999983E-003 s u(t)=  -1.0000000000000000E+000
 t=  4.3999999999999985E-003 s u(t)=  -1.0000000000000000E+000
 t=  4.4999999999999988E-003 s u(t)=  -1.0000000000000000E+000
 t=  4.5999999999999991E-003 s u(t)=  -1.0000000000000000E+000
 t=  4.6999999999999993E-003 s u(t)=  -1.0000000000000000E+000
 t=  4.7999999999999996E-003 s u(t)=  -1.0000000000000000E+000
 t=  4.8999999999999998E-003 s u(t)=  -1.0000000000000000E+000
 t=  5.0000000000000001E-003 s u(t)=  -1.0000000000000000E+000
 t=  5.1000000000000004E-003 s u(t)=  -1.0000000000000000E+000
 t=  5.2000000000000006E-003 s u(t)=  -1.0000000000000000E+000
 t=  5.3000000000000009E-003 s u(t)=  -1.0000000000000000E+000
 t=  5.4000000000000012E-003 s u(t)=  -1.0000000000000000E+000
 t=  5.5000000000000014E-003 s u(t)=  -1.0000000000000000E+000
 t=  5.6000000000000017E-003 s u(t)=  -1.0000000000000000E+000
 t=  5.7000000000000019E-003 s u(t)=  -1.0000000000000000E+000
 t=  5.8000000000000022E-003 s u(t)=  -1.0000000000000000E+000
 t=  5.9000000000000025E-003 s u(t)=  -1.0000000000000000E+000
 t=  6.0000000000000027E-003 s u(t)=  -1.0000000000000000E+000
 t=  6.1000000000000030E-003 s u(t)=  -1.0000000000000000E+000
 t=  6.2000000000000033E-003 s u(t)=  -1.0000000000000000E+000
 t=  6.3000000000000035E-003 s u(t)=  -1.0000000000000000E+000
 t=  6.4000000000000038E-003 s u(t)=  -1.0000000000000000E+000
 t=  6.5000000000000040E-003 s u(t)=  -1.0000000000000000E+000
 t=  6.6000000000000043E-003 s u(t)=  -1.0000000000000000E+000
 t=  6.7000000000000046E-003 s u(t)=  -1.0000000000000000E+000
 t=  6.8000000000000048E-003 s u(t)=  -1.0000000000000000E+000
 t=  6.9000000000000051E-003 s u(t)=  -1.0000000000000000E+000
 t=  7.0000000000000053E-003 s u(t)=  -1.0000000000000000E+000
 t=  7.1000000000000056E-003 s u(t)=  -1.0000000000000000E+000
 t=  7.2000000000000059E-003 s u(t)=  -1.0000000000000000E+000
 t=  7.3000000000000061E-003 s u(t)=  -1.0000000000000000E+000
 t=  7.4000000000000064E-003 s u(t)=  -1.0000000000000000E+000
 t=  7.5000000000000067E-003 s u(t)=  -1.0000000000000000E+000
 t=  7.6000000000000069E-003 s u(t)=  -1.0000000000000000E+000
 t=  7.7000000000000072E-003 s u(t)=  -1.0000000000000000E+000
 t=  7.8000000000000074E-003 s u(t)=  -1.0000000000000000E+000
 t=  7.9000000000000077E-003 s u(t)=  -1.0000000000000000E+000
 t=  8.0000000000000071E-003 s u(t)=  -1.0000000000000000E+000
 t=  8.1000000000000065E-003 s u(t)=  -1.0000000000000000E+000
 t=  8.2000000000000059E-003 s u(t)=  -1.0000000000000000E+000
 t=  8.3000000000000053E-003 s u(t)=  -1.0000000000000000E+000
 t=  8.4000000000000047E-003 s u(t)=  -1.0000000000000000E+000
 t=  8.5000000000000041E-003 s u(t)=  -1.0000000000000000E+000
 t=  8.6000000000000035E-003 s u(t)=  -1.0000000000000000E+000
 t=  8.7000000000000029E-003 s u(t)=  -1.0000000000000000E+000
 t=  8.8000000000000023E-003 s u(t)=  -1.0000000000000000E+000
 t=  8.9000000000000017E-003 s u(t)=  -1.0000000000000000E+000
 t=  9.0000000000000011E-003 s u(t)=  -1.0000000000000000E+000
 t=  9.1000000000000004E-003 s u(t)=  -1.0000000000000000E+000
 t=  9.1999999999999998E-003 s u(t)=  -1.0000000000000000E+000
 t=  9.2999999999999992E-003 s u(t)=  -1.0000000000000000E+000
 t=  9.3999999999999986E-003 s u(t)=  -1.0000000000000000E+000
 t=  9.4999999999999980E-003 s u(t)=  -1.0000000000000000E+000
 t=  9.5999999999999974E-003 s u(t)=  -1.0000000000000000E+000
 t=  9.6999999999999968E-003 s u(t)=  -1.0000000000000000E+000
 t=  9.7999999999999962E-003 s u(t)=  -1.0000000000000000E+000
 t=  9.8999999999999956E-003 s u(t)=  -1.0000000000000000E+000
 t=  9.9999999999999950E-003 s u(t)=  -1.0000000000000000E+000

Добавлено через 6 минут
Правильно ли считает в точке t= 1.0000000000000002E-003 s u(t)= 9.9999999999997835E-001 (как повысить универсальность программы, организовывать ввод сигналов, с использованием матпарсера, учесть границы предыдущего и следующего сигналов (учитывая принцип суперпозиции для линейных цепей))? Пример переделан под прямоугольные импульсы поочередно положительной и отрицательной полярности .
Подходят ли формулы под условия ( >=,<= , >,< ) ,кажется под новый сигнал после первого не включать ?

Pascal
1
2
3
4
5
  m:=0;
         if(t>=0   ) and   (t<= Signal[0].tend) then m:=0;
         if(t>= Signal[1].tbegin) and   (t<= Signal[1].tend) then m:=1;
         if(t>= Signal[2].tbegin) and   (t<= Signal[2].tend) then m:=2;
         if(t>= Signal[3].tbegin) and   (t<= Signal[3].tend) then m:=3;

Добавлено через 22 секунды
Или переставить границы?

USERPC01 · 02.09.2018, 04:47 **[ТС]**

Pascal
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
program paramtrapint;
 uses    Crt, Graph ,SysUtils ;
 
type
  //double=real;
  func = function(  x: double ): double;
   Signalarray=record
        Uinp : func;
             tbegin : double;
             tend: double;
     end;
 
 
  funcparamdint  =  function  ( tau: double ; t : double; eps  : double  ; const  Uin  : func  ; const h : func  ) : double;
 
 
 
 
 function  Derivative ( t : double; eps : double ; const   uinp : func     )   : double ;
begin
    Result:=  ( uinp(t+eps) -  uinp(t ) )/eps;
 
end;
 
 
 
function   deriveU  ( tau : double; eps1 : double ; const  Uin   :   func     )  : double;
 
begin
 Result:= Derivative (tau, 1e-8,    Uin     );
end;
 
 
 
function  f_tointegr( tau: double ; t : double; eps  : double  ;  const Uin  : func   ; const  h : func   ) : double;
 const
   eps1=1e-8;
begin
 Result:= deriveU(tau ,eps,  Uin    )*h(t-tau);
end;
 
 
 
 
 
 
 
 
 function trapparam (  a, b,t, eps, eps2: double;  const  intfunc: funcparamdint  ;  const uinp : func   ;  const h  : func   ): double;
var
  h1,s,s0,s1,sn ,fun1,fun2,fun3  : double;
  i,n: integer;
begin
 
s:=1;
sn:=101;
n:=4;
   fun1:=intfunc(a,t ,eps2  ,  uinp  ,  h     );
   fun2:=intfunc(b,t ,eps2  ,  uinp  ,  h    );
s0:=(fun1+fun2)/2;
 
s1:=intfunc ( ((a+b)/2),t ,eps2   ,   uinp  ,  h    );
while (Abs(s-sn)>eps) do
begin
 sn:=s;
 h1:=(b-a)/n;
 
    i:=0 ;
    while  i<(n/2)   do
    begin
       fun3:= intfunc( ( a+(2*i+1)*h1 ) , t  ,eps2 ,  uinp  , h   ) ;
    s1:=s1+fun3;
    i:=i+1;
    end;
 s:=h1*(s0+s1);
n:=n*2;
 end ;
Result:= s;
end;
     { **************************************** }
 
 
 
 
 
 
 
 
 
      function h( t :double) : double;
  var
    R,C,k,Tau : double;
 begin
 
   R:=1000;
   C:=0.1e-6;
   k:=1;   //0.25
   Tau:=R*C;
   if t<0 then Result:=0 else
    Result:=  k*(1-Exp(-t/Tau))   ;
 
 end;
 
 
 
    function Uinp1( t: double ) : double ;
   begin
 
      Result:=1 ;  //t 0...0.1 s
   end;
 
     function Uinp2( t: double ) : double ;
   begin
 
      Result:=0;   // 0.1 s...3 s
   end;
 
 
       function Uinp3( t: double ) : double ;
   begin
 
      Result:=-1; // 3s..3.2 s
   end;
 
     function Uinp4( t: double ) : double ;
   begin
 
      Result:=0;   // 3.2s..6 s
   end;
 
   function AssignSignal(const   Ui  : func; tbeg :double; tend: double ): Signalarray;
   var
     Uinput :Signalarray ;
   begin
 
              Uinput.Uinp:=Ui;
              Uinput.tbegin:=tbeg;
              Uinput.tend:=tend;
 
      Result:=Uinput;
   end;
 
 
 
  function getht( const h: func ; t : double)  : double;
begin
       Result:=h(t);
end;
 
 
 
 
 function DuhamelInt(   t: double ;  eps , eps2 : double  ;   const h:func ): double;
 
 
 
var
   dUinp_Ht  :  double;
   y: double;
   te  ,tp : double;
   Signal: array of Signalarray ;
   n,m ,i : integer;
 
 
begin
              //U(t):=U(0)*h(t)+  integr(0;t;derivU(tau)*h(t-tau); dtau)
//U(t):=U(0)*h(t)+  integr(0;t;derivU(tau)*h(t-tau); dtau)
 {
    u
 
 
    t=[0 ,t1]
   Y1(t) =U1(0)*h(t)+  integr(0;t;derivU(tau)*h(t-tau); dtau)  //0..t1
     Отклик на остальных интервалах вычисляется по формулам, вытекающих из принципа суперпозиции:
   Y2(t)=U1(0)*h(t)+integr(0;t1;derivU1(tau)*h(t-tau); dtau)+
   +(U2(t1)-U1(t1))*(h(t-t1))  +integr(t1; t ; derivU2(tau)*h(t-tau); dtau)+
 
   Y3= U1(0)*h(t)  +  integr(0;t1;derivU1(tau)*h(t-tau); dtau) +
      +(U2(t1)-U1(t1))*(h(t-t1)) +integr(t1; t2 ; derivU2(tau)*h(t-tau); dtau)+
      +(U3(t2)-U2(t2))*(h(t-t2)) +integr(t2; t  ; derivU3(tau)*h(t-tau); dtau)
 
}
        //for dummy  example , if U(t):=U(0)*h(t)+  integr(0;t;derivU(tau)*h(t-tau); dtau)
 
           n:=4;
           //n:=4
 
           SetLength(Signal,n);   //0,1,2- >U1,U2,U3
 
           Signal[0]:=AssignSignal(  @Uinp1  ,0, 1e-3  ) ;   //[0,t1]
           Signal[1]:=AssignSignal(  @Uinp2  ,1e-3,2e-3   ) ;   //[t1,t2]
           Signal[2]:=AssignSignal(  @Uinp3  ,2e-3,3e-3    ) ;   //[t2,t3]
           Signal[3]:=AssignSignal(  @Uinp4  ,3e-3,4e-3   ) ;   //[t3,t4]
 
 
 
 
 
          m:=0;
        // if(t>=0   ) and   (t<= Signal[0].tend) then m:=0;
       //  if(t>= Signal[1].tbegin) and   (t<= Signal[1].tend) then m:=1;
        // if(t>= Signal[2].tbegin) and   (t<= Signal[2].tend) then m:=2;
        // if(t>= Signal[3].tbegin) and   (t<= Signal[3].tend) then m:=3;
               if(t>=0   ) and   (t<= Signal[0].tend) then m:=0;
               if(t>=  Signal[1].tbegin) and   (t<= Signal[1].tend) then m:=1;
               if(t>= Signal[2].tbegin) and   (t<= Signal[2].tend) then m:=2;
               if(t>=   Signal[3].tbegin) and   (t<= Signal[3].tend) then m:=3;
               if  (t> Signal[3].tend)  then  begin writeln('end of signal'); Result:=0; exit; end;
 
 
 
 
               y:=0;
           for i:=0 to m do
            begin
 
 
               if (i=0) then
                 begin
 
 
 
                y :=y+Signal[0].Uinp(0)*getht(h,t)  +trapparam(0,t,t,eps,eps2, @f_tointegr,Signal[0].Uinp,  h    );      //U1
                 end
 
                  else
                    begin
                          tp:=Signal[i].tbegin  ;  //t1,t2 ;
                          te:=Signal[i].tend  ;    //t2,t3 ;
                          if (i=m)  then   te:=t;
                           dUinp_ht:=( Signal[i].Uinp(tp)-Signal[i-1].Uinp (tp) )*getht(h,t-tp);   //U2(t1)-U1(t1) ; U3(t2)-U2(t2) ;   U4(t3)-U3(t3)
                           y :=y+dUinp_ht+trapparam(tp,te, t,eps,eps2, @f_tointegr, Signal[i].Uinp , h  ); //U2,U3,U4
                     end;
            end;
                      Result:=y;
              // for Uinp=U0*1(t)   U(t)=U0 *(1-exp(-t/tau))
 
          // Result:=dUinp_Ht+trapparam(0,t, t,eps,eps2, @f_tointegr  ,  Uinp    ,  h   );
 
 
end;
 
 
 
 
 
 
 
 
 
 
 
 
procedure  main  ;
 var
      fp : TextFile;
      time,Tend,step,s: double;
begin
 
//writeln (' s(t)= trap(0,t,eps,f);   f(tau)=p(tau)*h(t-tau) ');
writeln ('  Input Tend, ms (6) ');
readln( Tend);
 Tend:=Tend*0.001 ;
writeln (' Input step ,ms,(0.1) ');
readln( step);
 
  step:=step*0.001;
 
     assignfile(fp, 'Integral.txt');
      ReWrite(fp);
 
 
writeln (fp,'  Tbeg=0 , Tend=', Tend);
writeln (fp,'   step=',step,'s ' );
//writeln (fp,'   p(tau)=-20000*exp(-200*tau); ');
//writeln (fp,'   h(t)=0.005-0.0025*exp(-100*t);  ');
 
  time:=0 ;
 
while    (time<=Tend) do
begin
 
s:=DuhamelInt(time, 1e-6,1e-8,   @h   );
 
 
 
 
writeln(' t= ',time,' s u(t)=  '  ,s) ;
writeln (fp,' t=',time, 's u(t)=  ' ,s);
time:=time+step  ;
end;
  CloseFile( fp);
readln;
 
end;
 
 
 
 
 
 
 
 
       {
 
 
 
 
 
  function trap (  a, b, eps: double; const intfunc: func  ): Double;
var
  h1,s,s0,s1,sn   : double;
  i,n: integer;
begin
 
s:=1;
sn:=101;
n:=4;
 
s0:=(intfunc(a )+intfunc(b ))/2;
 
s1:=intfunc ( ((a+b)/2)  );
while (Abs(s-sn)>eps) do
begin
 sn:=s;
 h1:=(b-a)/n;
 
    i:=0 ;
    while  i<(n/2)   do
    begin
 
    s1:=s1+intfunc( a+(2*i+1)*h1 ) ;  ;
    i:=i+1;
    end;
 s:=h1*(s0+s1);
n:=n*2;
 end ;
Result:= s;
end;
 
 
 
 
 
 
 
 
 
 function xsquare( x: double   ): double ;
begin
 
Result:=x*x ;
end;
 
 procedure GetIntegral;
 const
  eps=1e-8;
 var y: real;
 begin
    y:= trap(1,5, eps,@xsquare);
    writeln('y=',y);
    readln;
 end;
 
 
       }
 
 
 
 
 
 
 
 
            Function fgraph(x:real): real;
begin
Result:= DuhamelInt(x, 1e-6,1e-8,   @h   );
//fgraph:=100*sin(x) //сюда вписывать функцию графика
end;
 
 
 
  procedure Postroenie;  // строит график функции
var
x1,x2:real;    // границы изменения аргумента функции
y1,y2:real;    // границы изменения значения функции
x:real;        // аргумент функции
y:real;        // значение функции в точке x
step:real;     // приращение аргумента
l,b:integer;   // левый нижний угол области вывода графика
w,h:integer;   // ширина и высота области вывода графика
mx,my:real;    // масштаб по осям X и Y
x0,y0:integer; // точка - начало координат
n :integer;
str: string ;
gd,gm :smallint;
ymax ,ymin : real ;
begin
// область вывода графика
l:=100;
b:=600;
h:=300;          // высота
w:=1024;         // ширина
 
x1:=0;     // нижняя граница диапазона аргумента
x2:=0.0035;    // верхняя граница диапазона аргумента
step:=1e-8;  // шаг аргумента
 
 
 
 // вычислим масштаб
my:=100; // h/1e-12+abs(y2-y1);  // масштаб по оси Y
mx:=2*h/abs(x2-x1);  // масштаб по оси X
 
// оси
  x0:=l;
  y0:=b;
 
  Writeln('Initialising Graphics   .');
   gm:=0;
   gd:=0;
 
   initgraph(gd,gm,'');
 
// оси
 
  line(l,b+h,l ,b-h);
  line(l,b,l+w,b);
 
 
 
 
 
 
outtextXY(l+5,b-h,'Y');
outtextXY(l+w-10,b+10 ,'t');
 
// построение графика
x:=x1;
setcolor(12);
n:=0;
    ymax:=0;
    ymin:=0;
while ( x<x2) do
 begin
 
 y:=fgraph(x);
 
           if  ymax < y   then  ymax:=y;
           if  ymin > y      then  ymin:=y;
      PutPixel ((l+ Round(x*mx )) , (b- Round(y*my))   ,11)  ;
 
 
   // putpixel( x0+ Round(x*mx),b-Round(y*my),0);
       if (n=50000) then
        begin
            str:=FloatToStr (Round(x*100000)/100000)  ;
          outtextXY((l+ Round(x*mx )),b+50, str);
             str:=FloatToStr (Round(ymax*100000)/100000)  ;
            outtextXY (l-70,(b- Round(ymax*my)), str);
            str:=FloatToStr (Round(ymin*100000)/100000)  ;
            outtextXY (l-70,(b- Round(ymin*my)), str);
          n:=0;
          end;
 n:=n+1;
 
 x:=x + step;
end;
 
      readln;
     CloseGraph;
     Writeln('Graphics was Closed.');
end;
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
begin
  writeln('trap  test');
 // GetIntegral;
  writeln('trapparam test');
 // main;
 
  readln;
 
   Postroenie;
 
end.

Добавлено через 14 минут

Pascal
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
program paramtrapint;
 uses    Crt, Graph ,SysUtils ;
 
type
  //double=real;
  func = function(  x: double ): double;
   Signalarray=record
        Uinp : func;
             tbegin : double;
             tend: double;
     end;
 
 
  funcparamdint  =  function  ( tau: double ; t : double; eps  : double  ; const  Uin  : func  ; const h : func  ) : double;
 
 
 
 
 function  Derivative ( t : double; eps : double ; const   uinp : func     )   : double ;
begin
    Result:=  ( uinp(t+eps) -  uinp(t ) )/eps;
 
end;
 
 
 
function   deriveU  ( tau : double; eps1 : double ; const  Uin   :   func     )  : double;
 
begin
 Result:= Derivative (tau, 1e-8,    Uin     );
end;
 
 
 
function  f_tointegr( tau: double ; t : double; eps  : double  ;  const Uin  : func   ; const  h : func   ) : double;
 const
   eps1=1e-8;
begin
 Result:= deriveU(tau ,eps,  Uin    )*h(t-tau);
end;
 
 
 
 
 
 
 
 
 function trapparam (  a, b,t, eps, eps2: double;  const  intfunc: funcparamdint  ;  const uinp : func   ;  const h  : func   ): double;
var
  h1,s,s0,s1,sn ,fun1,fun2,fun3  : double;
  i,n: integer;
begin
 
s:=1;
sn:=101;
n:=4;
   fun1:=intfunc(a,t ,eps2  ,  uinp  ,  h     );
   fun2:=intfunc(b,t ,eps2  ,  uinp  ,  h    );
s0:=(fun1+fun2)/2;
 
s1:=intfunc ( ((a+b)/2),t ,eps2   ,   uinp  ,  h    );
while (Abs(s-sn)>eps) do
begin
 sn:=s;
 h1:=(b-a)/n;
 
    i:=0 ;
    while  i<(n/2)   do
    begin
       fun3:= intfunc( ( a+(2*i+1)*h1 ) , t  ,eps2 ,  uinp  , h   ) ;
    s1:=s1+fun3;
    i:=i+1;
    end;
 s:=h1*(s0+s1);
n:=n*2;
 end ;
Result:= s;
end;
     { **************************************** }
 
 
 
 
 
 
 
 
 
      function h ( t :double) : double;
  var
    R,C,k,Tau : double;
 begin
 
   R:=1000;
   C:=0.15e-6;
   k:=1;   //0.25
   Tau:=R*C;
   if t<0 then Result:=0 else
    Result:=  k*( Exp(-t/Tau))   ;
 
 end;
 
       function h2( t :double) : double;
  var
    R,C,k,Tau : double;
 begin
 
   R:=1000;
   C:=0.15e-6;
   k:=1;   //0.25
   Tau:=R*C;
   if t<0 then Result:=0 else
    Result:=  k*(1-Exp(-t/Tau))   ;
 
 end;
 
    function Uinp1( t: double ) : double ;
   begin
 
      Result:=1 ;  //t 0...0.1 s
   end;
 
     function Uinp2( t: double ) : double ;
   begin
 
      Result:=0;   // 0.1 s...3 s
   end;
 
 
       function Uinp3( t: double ) : double ;
   begin
 
      Result:=-1; // 3s..3.2 s
   end;
 
     function Uinp4( t: double ) : double ;
   begin
 
      Result:=0;   // 3.2s..6 s
   end;
 
   function AssignSignal(const   Ui  : func; tbeg :double; tend: double ): Signalarray;
   var
     Uinput :Signalarray ;
   begin
 
              Uinput.Uinp:=Ui;
              Uinput.tbegin:=tbeg;
              Uinput.tend:=tend;
 
      Result:=Uinput;
   end;
 
 
 
  function getht( const h: func ; t : double)  : double;
begin
       Result:=h(t);
end;
 
 function GetSignal(t: double  ):double;
 var
   Signal1: array of Signalarray ;
   m:integer;
 begin
      SetLength(Signal1,4);   //0,1,2- >U1,U2,U3
           m:=0;
           Signal1[0]:=AssignSignal(  @Uinp1  ,0, 1e-3  ) ;   //[0,t1]
           Signal1[1]:=AssignSignal(  @Uinp2  ,1e-3,2e-3   ) ;   //[t1,t2]
           Signal1[2]:=AssignSignal(  @Uinp3  ,2e-3,3e-3    ) ;   //[t2,t3]
           Signal1[3]:=AssignSignal(  @Uinp4  ,3e-3,4e-3   ) ;   //[t3,t4]
            if(t>=0   ) and   (t<= Signal1[0].tend) then m:=0;
               if(t>=  Signal1[1].tbegin) and   (t<= Signal1[1].tend) then m:=1;
               if(t>=  Signal1[2].tbegin) and   (t<= Signal1[2].tend) then m:=2;
               if(t>=    Signal1[3].tbegin) and   (t<= Signal1[3].tend) then m:=3;
 
             Result:=  Signal1[m].Uinp(t);
 
 end;
 
 function DuhamelInt(   t: double ;  eps , eps2 : double  ;   const h:func ): double;
 
 
 
var
   dUinp_Ht  :  double;
   y: double;
   te  ,tp : double;
   Signal: array of Signalarray ;
   n,m ,i : integer;
 
 
begin
              //U(t):=U(0)*h(t)+  integr(0;t;derivU(tau)*h(t-tau); dtau)
//U(t):=U(0)*h(t)+  integr(0;t;derivU(tau)*h(t-tau); dtau)
 {
    u
 
 
    t=[0 ,t1]
   Y1(t) =U1(0)*h(t)+  integr(0;t;derivU(tau)*h(t-tau); dtau)  //0..t1
     Отклик на остальных интервалах вычисляется по формулам, вытекающих из принципа суперпозиции:
   Y2(t)=U1(0)*h(t)+integr(0;t1;derivU1(tau)*h(t-tau); dtau)+
   +(U2(t1)-U1(t1))*(h(t-t1))  +integr(t1; t ; derivU2(tau)*h(t-tau); dtau)+
 
   Y3= U1(0)*h(t)  +  integr(0;t1;derivU1(tau)*h(t-tau); dtau) +
      +(U2(t1)-U1(t1))*(h(t-t1)) +integr(t1; t2 ; derivU2(tau)*h(t-tau); dtau)+
      +(U3(t2)-U2(t2))*(h(t-t2)) +integr(t2; t  ; derivU3(tau)*h(t-tau); dtau)
 
}
        //for dummy  example , if U(t):=U(0)*h(t)+  integr(0;t;derivU(tau)*h(t-tau); dtau)
 
           n:=4;
           //n:=4
 
           SetLength(Signal,n);   //0,1,2- >U1,U2,U3
 
           Signal[0]:=AssignSignal(  @Uinp1  ,0, 1e-3  ) ;   //[0,t1]
           Signal[1]:=AssignSignal(  @Uinp2  ,1e-3,2e-3   ) ;   //[t1,t2]
           Signal[2]:=AssignSignal(  @Uinp3  ,2e-3,3e-3    ) ;   //[t2,t3]
           Signal[3]:=AssignSignal(  @Uinp4  ,3e-3,4e-3   ) ;   //[t3,t4]
 
 
 
 
 
          m:=0;
        // if(t>=0   ) and   (t<= Signal[0].tend) then m:=0;
       //  if(t>= Signal[1].tbegin) and   (t<= Signal[1].tend) then m:=1;
        // if(t>= Signal[2].tbegin) and   (t<= Signal[2].tend) then m:=2;
        // if(t>= Signal[3].tbegin) and   (t<= Signal[3].tend) then m:=3;
               if(t>=0   ) and   (t<= Signal[0].tend) then m:=0;
               if(t>=  Signal[1].tbegin) and   (t<= Signal[1].tend) then m:=1;
               if(t>=  Signal[2].tbegin) and   (t<= Signal[2].tend) then m:=2;
               if(t>=    Signal[3].tbegin) and   (t<= Signal[3].tend) then m:=3;
               if  (t> Signal[3].tend)  then  begin writeln('end of signal'); Result:=0; exit; end;
 
 
 
 
               y:=0;
           for i:=0 to m do
            begin
 
 
               if (i=0) then
                 begin
 
 
 
                y :=y+Signal[0].Uinp(0)*getht(h,t)  +trapparam(0,t,t,eps,eps2, @f_tointegr,Signal[0].Uinp,  h    );      //U1
                 end
 
                  else
                    begin
                          tp:=Signal[i].tbegin  ;  //t1,t2 ;
                          te:=Signal[i].tend  ;    //t2,t3 ;
                          if (i=m)  then   te:=t;
                           dUinp_ht:=( Signal[i].Uinp(tp)-Signal[i-1].Uinp (tp) )*getht(h,t-tp);   //U2(t1)-U1(t1) ; U3(t2)-U2(t2) ;   U4(t3)-U3(t3)
                           y :=y+dUinp_ht+trapparam(tp,te, t,eps,eps2, @f_tointegr, Signal[i].Uinp , h  ); //U2,U3,U4
                     end;
            end;
                      Result:=y;
              // for Uinp=U0*1(t)   U(t)=U0 *(1-exp(-t/tau))
 
          // Result:=dUinp_Ht+trapparam(0,t, t,eps,eps2, @f_tointegr  ,  Uinp    ,  h   );
 
 
end;
 
 
 
 
 
 
 
 
 
 
 
 
procedure  main  ;
 var
      fp : TextFile;
      time,Tend,step,s: double;
begin
 
//writeln (' s(t)= trap(0,t,eps,f);   f(tau)=p(tau)*h(t-tau) ');
writeln ('  Input Tend, ms (6) ');
readln( Tend);
 Tend:=Tend*0.001 ;
writeln (' Input step ,ms,(0.1) ');
readln( step);
 
  step:=step*0.001;
 
     assignfile(fp, 'Integral.txt');
      ReWrite(fp);
 
 
writeln (fp,'  Tbeg=0 , Tend=', Tend);
writeln (fp,'   step=',step,'s ' );
//writeln (fp,'   p(tau)=-20000*exp(-200*tau); ');
//writeln (fp,'   h(t)=0.005-0.0025*exp(-100*t);  ');
 
  time:=0 ;
 
while    (time<=Tend) do
begin
 
s:=DuhamelInt(time, 1e-6,1e-8,   @h   );
 
 
 
 
writeln(' t= ',time,' s u(t)=  '  ,s) ;
writeln (fp,' t=',time, 's u(t)=  ' ,s);
time:=time+step  ;
end;
  CloseFile( fp);
readln;
 
end;
 
 
 
 
 
 
 
 
       {
 
 
 
 
 
  function trap (  a, b, eps: double; const intfunc: func  ): Double;
var
  h1,s,s0,s1,sn   : double;
  i,n: integer;
begin
 
s:=1;
sn:=101;
n:=4;
 
s0:=(intfunc(a )+intfunc(b ))/2;
 
s1:=intfunc ( ((a+b)/2)  );
while (Abs(s-sn)>eps) do
begin
 sn:=s;
 h1:=(b-a)/n;
 
    i:=0 ;
    while  i<(n/2)   do
    begin
 
    s1:=s1+intfunc( a+(2*i+1)*h1 ) ;  ;
    i:=i+1;
    end;
 s:=h1*(s0+s1);
n:=n*2;
 end ;
Result:= s;
end;
 
 
 
 
 
 
 
 
 
 function xsquare( x: double   ): double ;
begin
 
Result:=x*x ;
end;
 
 procedure GetIntegral;
 const
  eps=1e-8;
 var y: real;
 begin
    y:= trap(1,5, eps,@xsquare);
    writeln('y=',y);
    readln;
 end;
 
 
       }
 
 
 
 
 
 
 
 
            Function fgraph(x:real): real;
begin
Result:= DuhamelInt(x, 1e-6,1e-8,   @h   );
//fgraph:=100*sin(x) //сюда вписывать функцию графика
end;
 
  
            Function fgraph1(x:real): real;
begin
Result:= DuhamelInt(x, 1e-6,1e-8,   @h2   );
//fgraph:=100*sin(x) //сюда вписывать функцию графика
end;
              Function fgraph2(x:real): real;
begin
Result:=GetSignal(x  ) ;
//fgraph:=100*sin(x) //сюда вписывать функцию графика
end;
  procedure Postroenie;  // строит график функции
var
x1,x2:real;    // границы изменения аргумента функции
y1,y2:real;    // границы изменения значения функции
x:real;        // аргумент функции
y:real;        // значение функции в точке x
step:real;     // приращение аргумента
l,b:integer;   // левый нижний угол области вывода графика
w,h:integer;   // ширина и высота области вывода графика
mx,my:real;    // масштаб по осям X и Y
x0,y0:integer; // точка - начало координат
n :integer;
str: string ;
gd,gm :smallint;
ymax ,ymin : real ;
begin
// область вывода графика
l:=100;
b:=600;
h:=300;          // высота
w:=1024;         // ширина
 
x1:=0;     // нижняя граница диапазона аргумента
x2:=0.0035;    // верхняя граница диапазона аргумента
step:=1e-8;  // шаг аргумента
 
 
 
 // вычислим масштаб
my:=100; // h/1e-12+abs(y2-y1);  // масштаб по оси Y
mx:=2*h/abs(x2-x1);  // масштаб по оси X
 
// оси
  x0:=l;
  y0:=b;
 
  Writeln('Initialising Graphics   .');
   gm:=0;
   gd:=0;
 
   initgraph(gd,gm,'');
 
// оси
 
  line(l,b+h,l ,b-h);
  line(l,b,l+w,b);
 
 
 
 
 
 
outtextXY(l+5,b-h,'Y');
outtextXY(l+w-10,b+10 ,'t');
 
// построение графика
x:=x1;
setcolor(12);
n:=0;
    ymax:=0;
    ymin:=0;
while ( x<x2) do
 begin
 
 y:=fgraph(x);
 
           if  ymax  < y   then  ymax:=y;
           if  ymin > y      then  ymin:=y;
      PutPixel ((l+ Round(x*mx )) , (b- Round(y*my))   ,11)  ;
   y:=fgraph1(x);
 
           if  ymax  < y   then  ymax:=y;
           if  ymin > y      then  ymin:=y;
 
            PutPixel ((l+ Round(x*mx )) , (b- Round(y*my))   ,12)  ;
   y:=fgraph2(x);
 
           if  ymax  < y   then  ymax:=y;
           if  ymin > y      then  ymin:=y;
      PutPixel ((l+ Round(x*mx )) , (b- Round(y*my))   ,10)  ;
 
 
 
   // putpixel( x0+ Round(x*mx),b-Round(y*my),0);
       if (n=50000) then
        begin
            str:=FloatToStr (Round(x*100000)/100000)  ;
          outtextXY((l+ Round(x*mx )),b+50, str);
             str:=FloatToStr (Round(ymax*100000)/100000)  ;
            outtextXY (l-70,(b- Round(ymax*my)), str);
            str:=FloatToStr (Round(ymin*100000)/100000)  ;
            outtextXY (l-70,(b- Round(ymin*my)), str);
          n:=0;
          end;
 n:=n+1;
 
 x:=x + step;
end;
 
      readln;
     CloseGraph;
     Writeln('Graphics was Closed.');
end;
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
begin
  writeln('trap  test');
 // GetIntegral;
  writeln('trapparam test');
 // main;
 
  readln;
 
   Postroenie;
 
end.

USERPC01 · 06.09.2018, 00:33 **[ТС]**

Есть ли подпрограмма аппроксимации фрагмента входного сигнала интерполяционным полиномом на основе экспоненциальных рядов (моделирующие процессы, аналогичные процессам в апериодических и колебательно-апериодических звеньях ) , экспоненциально-гармоническими и экспоненциально-степенными рядами ( в некоторых случаях достаточно аппроксимации степенными полиномами Ньютона с переключателем интервалов , для которых коэффициенты справедливы)?

Новые блоги и статьи Все статьи Все блоги /
1С: Контроль уникальности заводского номера Maks 23.03.2026 Алгоритм контроля уникальности заводского (или серийного) номера на примере документа выдачи шин для спецтехники с табличной частью. Данные берутся из регистра сведений, по которому настроено. . .	Хочу заставить корпорации вкладываться в здоровье сотрудников: делаю мат модель здравосохранения anaschu 22.03.2026 e7EYtONaj8Y Z4Tv2zpXVVo https:/ / github. com/ shumilovas/ med2. git	1С: Программный отбор элементов справочника по группе Maks 22.03.2026 Установка программного отбора элементов справочника "Номенклатура" из модуля формы документа. В качестве фильтра для отбора справочника служит группа номенклатуры. Отбор по наименованию группы. . .	Как я обхитрил таблицу Word Alexander-7 21.03.2026 Когда мигает курсор у внешнего края таблицы, и нам надо перейти на новую строку, а при нажатии Enter создается новый ряд таблицы с ячейками, то мы вместо нервных нажатий Энтеров мы пишем любые буквы. . .
Krabik - рыболовный бот для WoW 3.3.5a AmbA 21.03.2026 без регистрации и смс. Это не торговля, приложение не содержит рекламы. Выполняет свою непосредственную задачу - автоматизацию рыбалки в WoW - и ничего более. Однако если админы будут против -. . .	1С: Программный отбор элементов справочника по значению перечисления Maks 21.03.2026 Установка программного отбора элементов справочника "Сотрудники" из модуля формы документа. В качестве фильтра для отбора служит значение перечислений. / / Событие "НачалоВыбора" реквизита на форме. . .	Переходник USB-CAN-GPIO Eddy_Em 20.03.2026 Достаточно давно на работе возникла необходимость в переходнике CAN-USB с гальваноразвязкой, оный и был разработан. Однако, все меня терзала совесть, что аж 48-ногий МК используется так тупо: просто. . .	Оттенки серого Argus19 18.03.2026 Оттенки серого Нашёл в интернете 3 прекрасных модуля: Модуль класса открытия диалога открытия/ сохранения файла на Win32 API; Модуль класса быстрого перекодирования цветного изображения в оттенки. . .

USERPC01 7 / 7 / 0 Регистрация: 29.06.2018 Сообщений: 1,536
	06.09.2018, 00:33 [ТС]
	Есть ли подпрограмма аппроксимации фрагмента входного сигнала интерполяционным полиномом на основе экспоненциальных рядов (моделирующие процессы, аналогичные процессам в апериодических и колебательно-апериодических звеньях ) , экспоненциально-гармоническими и экспоненциально-степенными рядами ( в некоторых случаях достаточно аппроксимации степенными полиномами Ньютона с переключателем интервалов , для которых коэффициенты справедливы)? 0