Re: Integrating a Piecewise function four times but different results

*To*: mathgroup at smc.vnet.net*Subject*: [mg121791] Re: Integrating a Piecewise function four times but different results*From*: Peter Pein <petsie at dordos.net>*Date*: Sun, 2 Oct 2011 02:37:08 -0400 (EDT)*Delivered-to*: l-mathgroup@mail-archive0.wolfram.com*References*: <j66e9t$hch$1@smc.vnet.net>

Am 01.10.2011 09:08, schrieb Rob Y. H. Chai: > Hello everyone, > > > > I tried to integrate a piecewise function 4 times but got different results > for two supposedly identical cases. The integration was done in the Mathematica > version 8.0.1.0. I am interested in what I did wrong, so any comments will > be appreciated. My code is below. > > > > Sincerely > > > > Rob Y. H. Chai > > UC Davis > > > > > > (* Introduction p1[x] as piecewise funtion *) > > p1[x_] = Piecewise[{{2*wmax*x/L, x <= L/2}, {2*wmax*(1 - x/L), x > L/2}}]; > > > > > > (* p2[x] is a specialized case of p1[x] with wmax = 1 and L = 1 *) > > p2[x_] = p1[x] /. wmax -> 1 /. L -> 1; > > > > > > (* now integrate p1[x] 4 times and introduce integration constants *) > > soln1[x_] = Integrate[Integrate[Integrate[Integrate[p1[x], x], x], x], x] + > a*x^3 + b*x^2 + c*x + d ; > > > > (* solve for integration constants in the solution of p1[x] *) > > eq1a = soln1[0] == 0; > > eq1b = (D[soln1[x], {x, 2}] /. x -> 0) == 0; > > eq1c = soln1[L] == 0; > > eq1d = ( D[soln1[x], {x, 2}] /. x -> L) == 0; > > soln1A[x_] = soln1[x] /. Simplify[Solve[{eq1a, eq1b, eq1c, eq1d}, {a, b, c, > d}], L > 0][[1]]; > > (* check case 1 solution at x = L/2 *) > > soln1A[L/2] /. wmax -> 1 /. L -> 1 > > 17/1920 > > > > > > (* similarly integrate p2[x] 4 times and introduce integration constants *) > > soln2[x_] = Integrate[Integrate[Integrate[Integrate[p2[x], x], x], x], x] + > r*x^3 + s*x^2 + u*x + v; > > (* solve for integration constants in the solution of p2[x] *) > > eq2a = soln2[0] == 0; > > eq2b = (D[soln2[x], {x, 2}] /. x -> 0) == 0; > > eq2c = soln2[1] == 0; > > eq2d = ( D[soln2[x], {x, 2}] /. x -> 1) == 0; > > soln2A[x_] = soln2[x] /. Simplify[Solve[{eq2a, eq2b, eq2c, eq2d}, {r, s, u, > v}]][[1]]; > > soln2A[1/2] (* case 2 solution at x = 1/2 *) > > 1/120 > > > > Plot[{soln1A[x] /. wmax -> 1 /. L -> 1, soln2A[x]}, {x, 0, 1}] > > > > Question - why doesn't case 1 solution converge to case 2 solution when case > 2 is simply a special case of case 1? Case 1 has got a discontinuity at L/2 (already after the first integration). I would prefer to iterate (Integrate w.r.t. x from 0 to z, add a constant and replace z by x) four times: assum = {L > 0 < wmax, Element[z, Reals]}; {p0[x_], p2[x_]} = FoldList[FullSimplify[PiecewiseExpand[ c[#2] + Integrate[#1, {x, 0, z}, Assumptions -> assum] /. z -> x], assum] & , Piecewise[{{(2*wmax*x)/L, x <= L/2}, {2*wmax*(1 - x/L), x > L/2}}], Range[4]] [[{-1, -3}]]; will give us the candidate function p0 and its second derivative. the conditions are: eqns = Simplify[Flatten[(#1 /@ {0, L} & ) /@ {p0, p2}] == 0, assum]; so we get: In[4]:= pfunc[x_] = FullSimplify[p0[x] /. First[ Solve[eqns, Union[Cases[eqns, _c, Infinity]]]], assum] Out[4]= Piecewise[{ {(wmax*x*(5*L^2 - 4*x^2)^2)/(960*L), x == 0 || L >= 2*x}}, (wmax*(L - x)*(L^2 + 8*L*x - 4*x^2)^2)/(960*L)] and the lim_{x->L/2}{pfunc(x)} is the same for both directions: In[5]:= (Limit[pfunc[x], x -> L/2, Assumptions -> assum, Direction -> #1] & ) /@ {1, -1} Out[5]= {(L^4*wmax)/120, (L^4*wmax)/120} Hope that helps, Peter