MathGroup Archive 1997

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: NIntegrate where terms of integrand have unknown constant coefficients

  • To: mathgroup at smc.vnet.net
  • Subject: [mg8964] Re: [mg8882] NIntegrate where terms of integrand have unknown constant coefficients
  • From: carlos at mars.Colorado.EDU (Carlos A. Felippa)
  • Date: Mon, 6 Oct 1997 01:59:33 -0400
  • Organization: University of Colorado, Boulder
  • Sender: owner-wri-mathgroup at wolfram.com

In article <611pv4$hlp at smc.vnet.net> David Withoff <withoff at wolfram.com> writes:
>>  I'm trying to do something along the lines of NIntegrate[c[1] f[1][x] +
>> c[2] f[2][x] + ..., {x, 0, a}], where c[n_] is unknown, but the f[n_] are
>> defined so that NIntegrate[f[n][x], {x, 0, a}] would work. I'd like to be
>> able to do the numerical integration, and keep the coefficients, so I'd get
>> as an answer c[1] NIntegrate[f[1][x], {x, 0, a}] + c[2] NIntegrate[f[2][x],
>> {x, 0, a}] + ... with all the NIntegrate's evaluated. The constant
>> coefficients are all of the same form c[n] (actually, each term will have
>> two coefficients c[n1] c[n2]).
>> Is there some way of doing this by changing the definition of NIntegrate so
>> it will automatically acheive this, or do I have to do something more
>> complicated? I'm not at all sure here, as I've never tried adding to the
>> definitions of complicated objects like NIntegrate.
>>
>> Any help or suggestions would be very gratefully accepted!
>>
>> Thanks, Scott Morrison
>> scott at morrison.fl.net.au
>
>I would do this by defining your own function that performs the
>symbolic linearity operations before calling NIntegrate, rather
>than by redefining NIntegrate.  For example
>
>In[1]:= int[p_Plus, q_] := Map[int[#, q] &, p]
>
>In[2]:= int[(p:c[_]) f_, q_] := p NIntegrate[f, q]
>
>In[3]:= int[c[1] x + c[2] x^2, {x, 0, 1}]
>
>Out[3]= 0.5 c[1] + 0.333333 c[2]
>
>This strategy could be made considerably more elaborate to do
>almost anything that you might want.
>
>Dave Withoff
>Wolfram Research
>

There is in fact another solution: write your own NIntegrate.
Below is an example that uses Gauss rules of order 1, 2 and 3.  The
logic can be easily extended to do variable intervals. For 
many problems that is not required if the symbolic part 
of f(x) is globally defined and the function is well behaved.

NIntGauss1[fx_,x_,lower_,upper_,n_]:= Module[{fint,d,x1,x2,f},
   	d=(upper-lower)/n; fint=0; f=fx;
   	Do[ x1=lower+(i-1)*d; x2=lower+i*d;  
	    fint +=(N[f/.x->N[(x1+x2)/2]])*N[d],
	   {i,1,n}];
	Return[fint]];
NIntGauss2[fx_,x_,lower_,upper_,n_]:= Module[{fint,d,x1,x2,f,
   	g1=N[(3+Sqrt[3])/6],g2=N[(3-Sqrt[3])/6]},
   	d=(upper-lower)/n; fint=0; f=fx;
   	Do[ x1=lower+(i-1)*d; x2=lower+i*d;  
        fint +=(N[f/.x->N[x1*g1+x2*g2]]+
                N[f/.x->N[x1*g2+x2*g1]])*N[d/2],
	   {i,1,n}];
	Return[fint]];
NIntGauss3[fx_,x_,lower_,upper_,n_]:= Module[{fint,d,x1,x2,f,
   	g1=N[(1+Sqrt[3/5])/2],g2=N[1/2],g3=N[(1-Sqrt[3/5])/2]},
   	d=(upper-lower)/n; fint=0; f=fx;
   	Do[ x1=lower+(i-1)*d; x2=lower+i*d;  
	    fint +=(N[5*f/.x->N[x1*g1+x2*g3]]+
                N[8*f/.x->N[(x1+x2)/2]]+
                N[5*f/.x->N[x1*g3+x2*g1]])*N[d/18],
	   {i,1,n}];
	Return[fint]];	
NIntGauss[f_,x_,lower_,upper_,n_,rule_]:= Module[{},
    If [rule==1,Return[NIntGauss1[f,x,lower,upper,n]]];
    If [rule==2,Return[NIntGauss2[f,x,lower,upper,n]]];
    If [rule==3,Return[NIntGauss3[f,x,lower,upper,n]]];
    Print["NintGauss: bad rule argument"];Return[Null]];

ClearAll[f]; 
f=a0+a1*x+a2*x^2+a3*x^3+a4*x^4+a5*x^5;
n=1; (* 3-pt Gauss is exact for any n, up to 5th order polys*)
Print[Simplify[NIntGauss[f,x,0,4,n,3]]];
f=Sin[x]^2*Sin[y]^2/N[Pi];
Print[Simplify[NIntGauss[f,x,0,Pi,4,3]]//InputForm];
Print[NIntegrate[f,{x,0,Pi}]];  (* fails *)

4. a0 + 8. a1 + 21.3333 a2 + 64. a3 + 204.8 a4 + 682.667 a5
0.5*Sin[y]^2



  • Prev by Date: Re: NIntegrate where terms of integrand have unknown constant coefficients
  • Next by Date: Parsing polygon lists
  • Previous by thread: Re: NIntegrate where terms of integrand have unknown constant coefficients
  • Next by thread: Re: FindRoot