U
USERPC01
Guest
You can translate into C++ this lazazrus script
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.tbegin ; //t1,t2 ;
te:=Signal.tend ; //t2,t3 ;
if (i=m) then te:=t;
dUinp_ht:=( Signal.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.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.
How to obtain Duhamel's integral using C/C++
Continue reading...
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.tbegin ; //t1,t2 ;
te:=Signal.tend ; //t2,t3 ;
if (i=m) then te:=t;
dUinp_ht:=( Signal.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.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.
How to obtain Duhamel's integral using C/C++
Continue reading...