{NotebookFileName[],DateString[]}

style[
t_
]:=
Style[t,Bold,FontFamily->"Arial",FontSize->14]

$Assumptions={{a,b,x,y}∈Reals}

The factorial

(∫_0^∞)(x^a)(E^(-x))\[DifferentialD]x

Assuming[
{a>0},
Factorial[a]==(∫_0^∞)(x^a)(E^(-x))\[DifferentialD]x//FullSimplify
]

Assuming[
{a>0},
Factorial[a+1]==(a+1)Factorial[a]//FullSimplify
]

(3/2)!

(1/2)!

(-(1/2))!

Gravitation due to a spherical shell

$Assumptions=Union[$Assumptions,{R>0,x>R,-R≤z≤R}]

We have a spherical shell of total mass M and radius R. A test mass m lies on the positive z-axis outside the sphere at z=x (thus x>R.) We want to find the gravitational force experienced by m due to the sphere. By symmetry, the force will be along the z-axis. We will find the total force by adding up the force from circular surface elements of depth dz, with z running from -R to R.

We decompose the shell into circular elements of depth dz. Here is an attempt to picture what one of these elements looks like.

Block[{is=400,R=1,z0=0.7,dz=0.05,r},
r=(Sqrt[(R^2)-(z0^2)]);
Show[
Graphics3D[{
{Opacity[0.2],
Sphere[{0,0,0},R]
},
{
Thickness[Scaled[0.01]],Opacity[1],Black,
Arrow[{{0,0,z0},{0,r,z0+dz}}],
Arrow[{{0,0,0},{0,0,z0}}],
Text[style["z"],{0,0,z0/2},{2,0}],
Text[style["r=(Sqrt[(R^2)-(z^2)])"],{0,r/2,z0+dz},{-1.1,0}],
Text[style["radius R sphere"],{0,0,-R}],
Text[style["radius r height dz surface element at z"],{0,-r,z0},{-0.3,1.5}]
}
}
,ImageSize->400
,Boxed->False
],
ContourPlot3D[(x^2)+(y^2)+(z^2)==1,{x,-1,1},{y,-1,1},{z,z0,z0+dz}
,PlotRange->{{-1,1},{-1,1},{-1,1}}
]
]
]

The element at z has radius r=(Sqrt[(R^2)-(z^2)]). Its width dw (this is the tricky part) extends from (z,r) to (z+dz,r+dr). To find dw we need dr in terms of z and dz. Circles will do for figuring this out:

Block[{R=1,z0=0.5,dz=0.05,φ0,dφ,r,dr,p1,p2},
r=(Sqrt[(R^2)-(z0^2)]);
dr=(Sqrt[(R^2)-(((z0+dz))^2)])-r;
φ0=ArcSin[(z0/R)];
dφ=ArcSin[((z0+dz)/R)]-φ0;
p1={r,z0};
p2={r+dr,z0+dz};
Graphics[{
{
Thick,Black,Circle[{0,0},R,{φ0,φ0+dφ}],
Text[style["dw"],((p1+p2)/2),{1,0.5}]
},
{Pink,Circle[{0,0},R,{0,(π/2)}]},
Arrow[{{0,0},p1}],Arrow[{{0,0},p2}],
{
Line[{p1,{r,z0+dz}}],
Text[style["dz"],{r,z0+(dz/2)},{-1.2,0}],
Line[{{r,z0+dz},p2}],
Text[style["dr"],{r+(dr/2),z0+dz},{0,-1}],
}
}
,PlotRange->{{0,1},{0,1}}
]
]

Zooming in:

Block[{R=1,z0=0.5,dz=0.05,φ0,dφ,r,dr,p1,p2},
r=(Sqrt[(R^2)-(z0^2)]);
dr=(Sqrt[(R^2)-(((z0+dz))^2)])-r;
φ0=ArcSin[(z0/R)];
dφ=ArcSin[((z0+dz)/R)]-φ0;
p1={r,z0};
p2={r+dr,z0+dz};
Graphics[{
{
Thick,Black,Circle[{0,0},R,{φ0,φ0+dφ}],
Text[style["dw"],((p1+p2)/2),{1,0.5}]
},
{Pink,Circle[{0,0},R,{0,(π/2)}]},
Arrow[{{0,0},p1}],Arrow[{{0,0},p2}],
{
Line[{p1,{r,z0+dz}}],
Text[style["dz"],{r,z0+(dz/2)},{-1.2,0}],
Line[{{r,z0+dz},p2}],
Text[style["dr"],{r+(dr/2),z0+dz},{0,-1}],
}
}
,PlotRange->{{r-0.5dr,r+1.5dr},{z0-0.5dz,z0+1.5dz}}
]
]

From this diagram we see that dw=(Sqrt[(((dz))^2)+(((dr))^2)]).

Since (z^2)+(r^2)=(R^2), we get by implicit differentiation 2zdz+2r dr=0, or dr=-(z/r)dz. Plugging in,

dw=(Sqrt[(((dz))^2)+((z^2)/(r^2))(((dz))^2)])
=(Sqrt[1+((z^2)/(r^2))])dz
=(Sqrt[(((r^2)+(z^2))/(r^2))])dz
=(Sqrt[((R^2)/(r^2))])dz
=(R/r)dz

Now, the area of the surface element is its circumference 2π r times its width dw, dA=2π rdw=2π Rdz, which is conveniently independent of z. (As a check, (∫_(z=-R)^R)dA=(∫_(-R)^R)2π Rdz=4π (R^2), as it should be. The mass of this surface element is thus (M/(4π (R^2)))2π Rdz=M((dz)/(2R)) (which also integrates up to M, as it should).

Block[{R=1,x=3,z0=0.5,r,angle=0.3},
r=(Sqrt[(R^2)-(z0^2)]);
Graphics[{
{
Thick, Black,PointSize[Large],
Arrow[{{0,0},{0,z0}}],
Arrow[{{0,x},{0,z0}}],
Circle[{0,z0},{r,angle r}],
Arrow[{{0,z0},{r,z0}}],
Text[style["z"],{0,(z0/2)},{-1.3,0}],
Text[style["x-z"],{0,((x+z0)/2)},{-1.3,0}],
Text[style["r"],{(r/2),z0},{1,-1.1}],
Point[{{0,0}}],
Text[style["center of sphere"],{0,0},{-1.1,0}],
Point[{{0,x}}],
Text[style["test mass"],{0,x},{-1.2,0}],
Arrow[{{0,x},{r,z0}}],
Text[style["force distance"],{(r/2),((x+z0)/2)},{-1.2,-2}],
Text[style["\!\(\*SqrtBox[
RowBox[{
StyleBox[SuperscriptBox[
RowBox[{"(", 
RowBox[{"x", "-", "z"}], ")"}], "2"],
FontSlant->"Italic"], 
StyleBox["+",
FontSlant->"Italic"], SuperscriptBox[
StyleBox["r",
FontSlant->"Italic"], "2"]}]]\)"],{(r/2),((x+z0)/2)},{-1.2,0.5}],
}
}
,PlotRange->{{-R,2R},{0,x}}
,PlotRangePadding->Scaled[.05]
]
]

The force on the test mass at x is a little more complicated. If all the mass of the surface element was axial, it would simply be -((G M m)/((((x-z))^2)+(r^2))). (Note: the distance from the test mass to any little bit of the surface element is (Sqrt[(((x-z))^2)+(r^2)]), which is why (((x-z))^2)+(r^2) appears in the denominator of the inverse square force.) But the mass is all off-axis by a distance r. The axial component of the force is

(Subscript[F, z])=(Subscript[F, total])(((x-z)/(Sqrt[(((x-z))^2)+(r^2)])))					
=-(((G M m)/(2R)))(1/((((x-z))^2)+(r^2)))((x-z)/(Sqrt[(((x-z))^2)+(r^2)]))
=-((G M m)/(2R))((x-z)/((((((x-z))^2)+(r^2)))^(3/2)))
=-((G M m)/(2R))((x-z)/((((x^2)-2x z+(z^2)+(r^2)))^(3/2)))
=-((G M m)/(2R))((x-z)/((((x^2)-2x z+(R^2)))^(3/2)))

(∫_(-R)^R)((x-z)/((((x^2)-2x z+(R^2)))^(3/2)))\[DifferentialD]z//Simplify

∫((x-z)/((((x^2)-2x z+(R^2)))^(3/2)))\[DifferentialD]z//Simplify

∫(x/((((x^2)-2x z+(R^2)))^(3/2)))\[DifferentialD]z//Simplify

∫((-z)/((((x^2)-2x z+(R^2)))^(3/2)))\[DifferentialD]z

∫(1/(x(Sqrt[(x^2)-2x z+(R^2)])))\[DifferentialD]z

How do I do that second integral? By parts, maybe?

u=z,du=dz
dv=((dz)/((((x^2)-2x z+(R^2)))^(3/2))),v=(1/(x(Sqrt[(x^2)-2x z+(R^2)])))

I=(z/(x(Sqrt[(x^2)-2x z+(R^2)])))-∫((dz)/(x(Sqrt[(x^2)-2x z+(R^2)])))
=(z/(x(Sqrt[(x^2)-2x z+(R^2)])))+((Sqrt[(x^2)-2x z+(R^2)])/(x^2))
=((x z+(x^2)-2x z+(R^2))/((x^2)(Sqrt[(x^2)-2x z+(R^2)])))
=(((x^2)-x z+(R^2))/((x^2)(Sqrt[(x^2)-2x z+(R^2)])))

Great! And now I just have to subtract that from (1/(Sqrt[(x^2)-2x z+(R^2)])), giving

∫((x-z)/((((x^2)-2x z+(R^2)))^(3/2)))dz=((x^2)/((x^2)(Sqrt[(x^2)-2x z+(R^2)])))-(((x^2)-x z+(R^2))/((x^2)(Sqrt[(x^2)-2x z+(R^2)])))
=((x z-(R^2))/((x^2)(Sqrt[(x^2)-2x z+(R^2)])))

At z=R the radical in the denominator simplifies to (Sqrt[(x^2)-2x R+(R^2)])=(Sqrt[(((x-R))^2)])=x-R, and the numerator is x R-(R^2)=R(x-R), so it all simplifies to (R/(x^2)). At z=-R we end up with ((-x R-(R^2))/((x^2)(x+R)))=-(R/(x^2)), and we’re there.

((x z-(R^2))/((x^2)(Sqrt[(x^2)-2x z+(R^2)])))/.{{z->R},{z->-R}}

Could I just do the whole thing by parts, without separating the terms?

u=x-z,du=-dz
					dv=((dz)/((((x^2)-2x z+(R^2)))^(3/2))),v=(1/(x(Sqrt[(x^2)-2x z+(R^2)])))

I=((x-z)/(x(Sqrt[(x^2)-2x z+(R^2)])))-∫((-dz)/(x(Sqrt[(x^2)-2x z+(R^2)])))
=((x-z)/(x(Sqrt[(x^2)-2x z+(R^2)])))-((Sqrt[(x^2)-2x z+(R^2)])/(x^2))
=(((x^2)-x z-((x^2)-2x z+(R^2)))/((x^2)(Sqrt[(x^2)-2x z+(R^2)])))
=((x z-(R^2))/((x^2)(Sqrt[(x^2)-2x z+(R^2)])))

Yes, that works, and is simpler.

This even works for |x|<R and x<-R. When x<R, x-R<0, so at z=R we have

I=((x R-(R^2))/((x^2)(Sqrt[(x^2)-2x R+(R^2)])))=((R(x-R))/((x^2)(Sqrt[(((x-R))^2)])))=((R(x-R))/(-(x^2)(x-R)))=-(R/(x^2))

The minus sign comes because, for negative x-R, (Sqrt[(((x-R))^2)])=-(x-R). x+R is still positive in this range, so at the lower bound z=-R we have

I=((-x R-(R^2))/((x^2)(Sqrt[(x^2)+2x R+(R^2)])))=((-R(x+R))/((x^2)(x+R)))=-(R/(x^2))

So the integral evaluates to I(Subscript[|, -R]^R)=-(R/(x^2))-(-(R/(x^2)))=0.

For x<-R both x-R and x+R are negative, so we end up with the negative of the force we got for x>R, which is appropriate, since the force is now in the opposite direction, pointing up instead of down.

Integration by parts examples

∫x (E^x)\[DifferentialD]x

∫(E^x)Sin[x]\[DifferentialD]x

∫(x^n)(E^(-x))\[DifferentialD]x

∫ArcTan[x]\[DifferentialD]x

∫((Sin[x])^n)\[DifferentialD]x

∫x (E^(-(x^2)/2))\[DifferentialD]x

∫(x^2) (E^(-(x^2)/2))\[DifferentialD]x

(∫_(-∞)^∞)(x^2)(E^(-(x^2)/2))\[DifferentialD]x

I can probably get that by parts.

u=x,du=dx
dv=x (e^(-(x^2)/2))dx, v=-(e^(-(x^2)/2))

(∫_(-∞)^∞)(x^2)(e^(-(x^2)/2))dx=x (e^(-(x^2)/2))(Subscript[|, -∞]^∞)-(∫_(-∞)^∞)(-(e^(-(x^2)/2)))dx
=0+(∫_(-∞)^∞)(e^(-(x^2)/2))dx

$Assumptions=Union[$Assumptions,{σ>0,μ∈Reals}]

(∫_(-∞)^∞)(E^(-((((x-μ))^2)/(2(σ^2)))))\[DifferentialD]x

(∫_(-∞)^∞)(1/((Sqrt[2π])σ))(E^(-((((x-μ))^2)/(2(σ^2)))))\[DifferentialD]x

(∫_(-∞)^∞)x(1/((Sqrt[2π])σ))(E^(-((((x-μ))^2)/(2(σ^2)))))\[DifferentialD]x

If I sub u=x-μ, this becomes

(∫_(-∞)^∞)(u+μ)(e^(-((u^2)/(σ^2))))du=0+μ

How about subbing u=((x-μ)/σ)? Then x=σ u+μ, dx=σdu, and for the 0-2 moments I get

(1/((Sqrt[2π])σ))(∫_(-∞)^∞)(e^(-((u^2)/2)))σdu=(1/(Sqrt[2π]))(∫_(-∞)^∞)(e^(-((u^2)/2)))du=1
(1/((Sqrt[2π])σ))(∫_(-∞)^∞)(σ u+μ)(e^(-((u^2)/2)))σdu=(1/(Sqrt[2π]))(∫_(-∞)^∞)μ (e^(-((u^2)/2)))du=μ

(1/((Sqrt[2π])σ))(∫_(-∞)^∞)(((σ u+μ))^2)(e^(-((u^2)/2)))σdu=(1/(Sqrt[2π]))(∫_(-∞)^∞)(((σ u+μ))^2)(e^(-((u^2)/2)))du
=(1/(Sqrt[2π]))(∫_(-∞)^∞)((σ^2) (u^2)+2σ μ u+(μ^2))(e^(-((u^2)/2)))du
=(σ^2)+0+(μ^2)=(σ^2)+(μ^2)

(∫_(-∞)^∞)(x^2)(1/((Sqrt[2π])σ))(E^(-((((x-μ))^2)/(2(σ^2)))))\[DifferentialD]x

∫(x^2)Cos[x]\[DifferentialD]x

∫(x/(Sqrt[x+1]))\[DifferentialD]x

∫(x/(Sqrt[x+c]))\[DifferentialD]x

∫(x/(Sqrt[c x+1]))\[DifferentialD]x

D[(2/c)(Sqrt[c x+1]),x]

D[((2c x-4)/(3(c^2)))(Sqrt[c x+1]),x]==(x/(Sqrt[c x+1]))//Simplify

∫Cos[x](E^x)\[DifferentialD]x

Path integral of a parabola

∫(Sqrt[1+(((2a t+b))^2)])\[DifferentialD]t

∫(Sqrt[1+(((2a t))^2)])\[DifferentialD]t

((Sinh[t])^2)+1//FullSimplify

So i guess we do this by substituting t=(1/(2a))sinh(z). Then I end up needing to integrate

(Sqrt[1+(((2a t))^2)])/.t->((Sinh[z])/(2a))

(Sqrt[1+(((2a t))^2)])/.t->((Sinh[z])/(2a))//FullSimplify

D[Sinh[z],z]

So, I end up integrating (cosh^2)z dz. There’s probably a double identity for (cosh^2):

((Cosh[z])^2)//TrigToExp

Cosh[2z]//TrigToExp

(1/2)(Cosh[2z]+1)//TrigToExp

((Cosh[z])^2)==(1/2)(Cosh[2z]+1)//Simplify

Right. That’s where the two terms come from. Wonder if I could also get them from IBP in some way? Kind of looks like it might be possible: the (1/2)t(Sqrt[1+4(a^2)(t^2)]) term is a product.

OK, let’s look at the “base case”, ((x^2)/2).

∫(Sqrt[1+(t^2)])\[DifferentialD]t

D[(Sqrt[1+(t^2)]),t]

u=(Sqrt[1+(t^2)]), du=(t/(Sqrt[1+(t^2)]))dt
dv=dt,v=t
∫(Sqrt[1+(t^2)])dt=t(Sqrt[1+(t^2)])-∫((t^2)/(Sqrt[1+(t^2)]))dt

That doesn’t seem to have helped much. But now a substitution might work:

t=sinh(z), dt=cosh(t)dz

D[Sinh[z],z]

((Sinh[z])^2)+1//FullSimplify

∫(((sinh^2)(z))/(cosh(z)))cosh(z)dz=∫(sinh^2)(z)dz=(1/2)∫(cosh(2z)-1)dz

((Sinh[z])^2)==(1/2)(Cosh[2z]-1)//Simplify

So that will work, but the IBP didn’t actually simplify things at all.

D[t(Sqrt[1+(t^2)]),t]

∫((Sec[θ])^3)\[DifferentialD]θ

(Sqrt[1+(x^2)])/.x->Tan[θ]//FullSimplify

1+(t^2)/.t->(1/2)(u-(1/u))//Simplify

(1+(t^2)/.t->(1/2)(u-(1/u)))==((((1/2)(u+(1/u))))^2)//Simplify

Solve[t==(1/2)(u-(1/u)),u]

$Assumptions=Union[$Assumptions,{t∈Reals}]

t+(Sqrt[1+(t^2)])>0//FullSimplify

If I make the substitution t=(1/2)(u-(1/u)), dt=(1/2)(1+(1/(u^2)))du. Then I get

∫(Sqrt[1+(t^2)dt])=∫(1/2)(u+(1/u))(1/2)(1+(1/(u^2)))du=(1/4)∫((1/(u^3))+(2/u)+u)du

(u+(1/u))(1+(1/(u^2)))//Expand

(1/4)∫(u+(1/u))(1+(1/(u^2)))\[DifferentialD]u

((1/4)∫(u+(1/u))(1+(1/(u^2)))\[DifferentialD]u)/.u->t+(Sqrt[1+(t^2)])

((1/4)∫(u+(1/u))(1+(1/(u^2)))\[DifferentialD]u)/.u->t+(Sqrt[1+(t^2)])//FullSimplify

Yup, that works. So I can do this integral by substitution without using a hyperbolic function.

Radioactive decay

$Assumptions=Union[$Assumptions,{th>0}]

p1[
t_
]:=
(2^(-t/th))

∫p1[t]\[DifferentialD]t

(∫_0^∞)p1[t]\[DifferentialD]t

(∫_0^∞)t p1[t]\[DifferentialD]t

(∫_0^∞)(t^2) p1[t]\[DifferentialD]t

(∫_0^∞)(t^n) p1[t]\[DifferentialD]t