# How to obtain Duhamel's integral using C/C++ • ### 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 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 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