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

  • 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