# Alternative version of program for Duhamel's integral parsing using C++ with derivative parser for Uin(t)

• ### Question

• Alternative version of program  for Duhamel's integral parsing using C++  with derivative parser for Uin(t) . Some bugs are fixed.

#include <stdio.h> //cstdio.h

#include <math.h>//cmath.h
#include <conio.h>//

/*
ПХ ,1(t)
Uout(t)=integr(0;t; U'inp(tau) * h(t-tau);d tau   )

U'inp(tau)= U(tau)-U(tau0)/(tau-tau0)
ИХ, delta(t)
Uout(t)=integr(0;t; Uinp(tau) * k(t-tau);d tau   )
k(t)=dh(t)/dt
*/

double Heaviside(double t)
{

double y;
if (t<0) {y= 0;}
if (t>=0) { y= 1;}
//printf("\n ***  arg=%le Uinp(t)=%le",t,y);
//getch();
return y;
}

double delta(double t)
{
if (t=0){ return 1;} else {return 0;}
}

double function(double t)
{

return Heaviside(t) ;
}

double getdydt(double tau , double (*Uin )(double )  )
{
double delta= 1e-10;
//printf("\n tau=%le Uin ( tau )=%le\n tau+delta=%le  Uin(tau+delta)=%le  ", tau, Uin(tau), tau+delta, Uin(tau+delta)     );

double y=(Uin(tau+delta)-Uin(tau )) /delta;
return y;
}

double p(double tau, double (*Uin )(double  ) )
{
/*dU/dt, t=tau;*/
return  getdydt(tau,Uin);
}

double taurc;

double h(double tau)
{

double R=100;
double C=1E-6;
double k=1;
double Tau=R*C;
taurc=Tau;
return k*(1-exp(-tau/Tau));

}

double f(double tau, double t ,double du )
{
//p(tau)=dUinp(tau)/d tau  Duhamel, 1(t)
// p(tau)=Uinp(tau) , delta(t)
//double    Uinp (double tau);
//p(tau,Uinp );
return  du*h(t-tau);
}

double trap(double a,double b,double eps, double(*f)(double  , double  , double U ), double (*Uin)(double tau )    )
{

double tau;
double h1,s,s0,s1,sn, t, y1,y2,U;
int i,n;
s=1;
sn=101;
n=4;
tau=a;
t=b;
U=p(tau, Uin );
y1=f(tau,t, U);
tau=b;
U=p(tau, Uin );
y2=f(tau,t, U);

s0=(y1+y2)/2;
tau=((a+b)/2);
U=p(tau, Uin );
s1=f(tau,t,U);

while (fabs(s-sn)>eps){
sn=s;
h1=(b-a)/n;
for (i=0; i<n/2; i++)
{
tau=(a+(2*i+1)*h1);
U=p(tau, Uin );
s1+=f(tau,t, U );
}
s=h1*(s0+s1);
n*=2;
}
return s;
}

double DuhamelInt( double t)
{
double a,b,er,eps;

//f(double tau, double t)
double f(double, double, double );
double s;
double trap(double,double,double, double(*f) (double,double , double  ), double (*) (double ) );

eps=0.0000001;

return function(0)*h(t)*1 +trap(0,t,eps,f, function );

//function(0)*h(t)* 0 for testing only trap() parser
}

int main()
{
FILE *fp;
double  time,Tend,step,s;
printf ("\n s(t)= trap(0,t,eps,f); \n f(tau)=Uin'(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,Heaviside); \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)=  ; ");
fprintf (fp,"\n   h(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);
}
printf ("\n tau RC=%le s" ,taurc );

for(time=0; time<=Tend; time=time+step)

s=getdydt(time , function  );

printf ("\n t= %le s u(t)=%le, U'(t)=%le ",time,function(time),s);
}
printf ("\n tau RC=%le s" ,taurc );

fclose(fp);

getch() ;
}

Friday, November 25, 2016 11:39 AM

### All replies

• How to use  double(*f) (double,double , double (*) (double )  ) in

double trap(double,double,double, double(*f) (double,double , double  ), double (*) (double ) );

for  make it  double trap(double,double,double, double(*f) (double,double , double (*) (double )  ), double (*) (double ) );  for  use and input  fuction , for example, Heaviside(t) as function-prototype in integration suboutine for using it as method of class?

Friday, November 25, 2016 11:46 AM