locked
How to obtain Duhamel's integral using C/C++ RRS feed

  • Question

  • #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() ;
    }


    • Edited by USERPC01 Saturday, February 7, 2015 2:19 PM
    Saturday, February 7, 2015 2:19 PM

All replies

  •  How to  declare  in VBA  double trap(double a,double b,double eps, double(*f) (double , double )) as Function (   ByVal a As Double, ByVal b As Double ,  ByVal eps As Double,  ByRef   f(  ? ) As Double   )     

    witout using this method  as in MyFx()  (   ByRef   f(  ? ) As Double )    ?

    How to patch "Function Romberg(ByVal a As Double, ByVal b As Double, ByVal eps As Double)" to load name of user defined function fun(x1,x2,x3,xn)  , argument (dx)  and vector of parameters( x2,x3,xn if it is necessary )    ) ?

    'Option Explicit
     
    ' Function MyFx(ByVal sExpress As String, ByVal sVarName As String, ByVal X As Double) As Double
    'MyFx = Evaluate(Replace(sExpress, sVarName, "(" & CStr(X) & ")"))
    'End Function       
     
     
     
     
     
    Function fun(ByVal x As Double) As Double
    fun = x * x

    End Function

    Function Romberg(ByVal a As Double, ByVal b As Double, ByVal eps As Double) As Double
     Dim h1 As Double
     Dim s  As Double
     Dim s0 As Double
     Dim s1 As Double
     Dim sn As Double
     Dim t As Double
     Dim i As Integer
     Dim n As Integer
     Dim fa As Double
     Dim fb As Double
     Dim fhab  As Double
     Dim ytemp As Double
     Dim x As Double
     s = 1
     sn = 101
     n = 4
     x = a
     fa = f(x)
     x = b
     fb = f(x)
     x = (a + b) / 2
     fhab = f(x)
     
     s0 = (fa + fb) / 2
     s1 = fhab
     
     Do While Abs(s - sn) > eps
       sn = s
       h1 = (b - a) / n
        i = 0
        Do While (i < n / 2)
          x = (a + (2 * i + 1) * h1)
          ytemp = f(x)
          s1 = s1 + ytemp
        Loop
       
          s = h1 * (s0 + s1)
          n = n * 2
        i = i + 1
     Loop

    Romberg = s
    End Function

    Sub Doint()

    Dim y As Double
    Dim eps As Double
    eps = 0.1

    'y = fun(2)  ' MsgBox (y)

     y = Romberg(1, 5, eps, fun)
    MsgBox (y)

    End Sub



    • Edited by USERPC01 Thursday, November 12, 2015 4:24 PM
    Thursday, November 12, 2015 4:17 PM
  • See Namir Clement Shammas A new face of of Romberg Integration :

        

    Function MyFx(ByVal sExpress As String, ByVal sVarName As String, ByVal X As Double) As Double

    MyFx = Evaluate(Replace(sExpress, sVarName, "(" & CStr(X) & ")"))

    End Function

    Function RombergBasic(ByVal sExpress As String, ByVal sVarName As String, _

    ByVal A As Double, ByVal B As Double, ByVal Toler As Double) As Double

    ' NOTE: The constants ROW0 and COL0 are offset values used to determine the lower values

    ' for the row and column indices of matrix R(,). These offsets are helpful when translating

    ' the VBA function into other BASIC dialects or languages that do not support zero indices

    ' for the rows and columns of a matrix.

    Const ROW0 = 1

    Const COL0 = 1

    Dim I As Integer, J As Integer, MaxCols As Integer, M As Long

    Dim R() As Double, h As Double, X As Double, Sum As Double

    Dim Func As String

    MaxCols = CInt(Abs(Log(Toler) / Log(10)))

    ReDim R(1 + ROW0, MaxCols + COL0)

    h = B - A

    R(ROW0, COL0) = h / 2 * (MyFx(sExpress, sVarName, A) + _

    MyFx(sExpress, sVarName, B))

    For I = 1 To MaxCols

    h = h / 2

    Sum = 0

    For J = 1 To 2 ^ (I - 1)

    Sum = Sum + MyFx(sExpress, sVarName, A + (2 * J - 1) * h)

    Next J

    R(1 + ROW0, COL0) = R(ROW0, COL0) / 2 + h * Sum

    M = 1

    For J = 1 To I

    M = 4 * M

    R(1 + ROW0, J + COL0) = (M * R(1 + ROW0, J - 1 + COL0) - _

    R(0 + ROW0, J - 1 + COL0)) / (M - 1)

    Next J

    For J = 0 To I

    R(0 + ROW0, J + COL0) = R(1 + ROW0, J + COL0)

    Next J

    Next I

    RombergBasic = R(0 + ROW0, MaxCols + COL0)

    End Function

    Copyright © 2012, 2014 Namir Clement Shammas

     

    Thursday, November 12, 2015 4:27 PM