Numerical Computing: Advanced Techniques

Applications and tips

Andrew Moylan
Wolfram Research
training@wolfram.com

Numerical Integration: Quadrature

NIntegrate does numerical integration of functions (“numerical quadrature”):

NIntegrate[2(Sqrt[1-(y^2)]),{y,-1,1}]

Plot[(Sqrt[1-(y^2)]),{y,-1,1},Filling->Axis,AspectRatio->Automatic]

Numerical Integration: List Integration

For numerical integration of a list of points, interpolate the points.

y={{0,1.01},{0.2,1.03},{0.4,0.93},{0.6,0.85},{0.8,0.72},{1,0.57}};

f=Interpolation[y];

Then use NIntegrate:

NIntegrate[f[t],{t,0,1}]

Clear[f]

Numerical Integration: Path Integration

For path (contour) integration, parameterize the path.

path[t_]:=2+Exp[I t]

f[z_]:=((Sin[z])/((z-2)^3))

Path integral:

NIntegrate[f[path[t]] path'[t],{t,0,2π}]

2π I Residue[f[t],{t,2}]//N

Numerical Integration: Vector Calculus

Surface integrals, volume integrals, and vector calculus.

Parametric surface:

{x,y,z}={(2+Cos[v])Cos[u],(2+Cos[v])Sin[u],Sin[v](5/4-Cos[u])};

ParametricPlot3D[{x,y,z},{u,0,2π},{v,0,2π}]

Surface area integral:

NIntegrate[Norm[Cross[D[{x,y,z},u],D[{x,y,z},v]]],{u,0,2π},{v,0,2π}]

Compute enclosed volume using the divergence theorem (Gauss’ theorem):

NIntegrate[(({x,y,z})/3).Cross[D[{x,y,z},u],D[{x,y,z},v]],{u,0,2π},{v,0,2π}]

Clear[x,y,z]

Differential Equations

NDSolve solves ordinary and partial differential equations.

NDSolve[{x''[t]==-x[t],x'[0]==x[0]==1},x,{t,0,10}]

The solution is a list of rules:

solution=x[t]/.First[%]

Plot[solution,{t,0,10}]

Differential Equations: Special Kinds of ODEs

NDSolve handles delay and differential-algebraic equations.

Delay differential equation:

{x[t],x'[t]}/.First[NDSolve[{x'[t]==x[t](x[t-Pi]-x'[t-1]),x[t/;t≤0]==Cos[t]},x,{t,0,5}]]

Plot[%,{t,0,5}]

Differential Equations: PDEs

Partial differential equation (Schrödinger equation):

equations={I D[ψ[x,t],t]==-(1/2)D[ψ[x,t],{x,2}]+250Sin[x/2]^10ψ[x,t],ψ[x,0]==Exp[I 20 x]Exp[-10x^2],ψ[-π,t]==ψ[π,t]};

sol=ψ/.First[NDSolve[equations,ψ,{x,-π,π},{t,0,0.4}]]

It’s a particle rolling back and forth in a periodic well:

Animate[Plot[{Sin[(x/2)]^10+Abs[sol[x,t]]^2,Sin[(x/2)]^10},{x,-π,π},PlotRange->{0,3/2},Filling->1->{2}],{t,0,0.4},DefaultDuration->10]

Differential Equations: Indefinite Integral

Use NDSolve to get a numerical indefinite integral.

Indefinite integral as a differential equation:

f=y/.First[NDSolve[{y'[x]==Sin[Exp[x]],y[0]==0},y,{x,0,3}]]

Plot the integral and integrand:

Animate[Plot[{f[x],Sin[Exp[x]]},{x,0,t},Filling->{2->Axis},PlotStyle->{Thick,None},PlotRange->{{0,3},{-1,1}}],{t,0,3}]

Equation Solving

FindRoot searches for a numerical solution of one or more equations.

FindRoot[Cos[Cos[x]]==x,{x,0}]

Equation Solving

You can use FindRoot to set up an inverse function.

Given y(x), find the inverse function inv(t):

y[x_]:=x^x+x

inv[t_]/;NumericQ[t]:=x/.FindRoot[y[x]==t,{x,1}]

Use it like any other function:

Plot[inv[t],{t,1,10}]

Optimization

Local and global optimization with constraints.

FindMinimum looks for a nearby local minimum:

FindMinimum[Exp[-x^2] Cos[5Exp[x]],{x,1}]

NMinimize works harder to find the global minimum:

NMinimize[Exp[-x^2] Cos[5Exp[x]],x]

Optimization

Multidimensional optimization with constraints:

NMinimize[{Sin[x]Cos[2y+1],x^4+y^2<3},{x,y}]

Linear Algebra: Linear Solving

LinearSolve is the workhorse in many numerical algorithms.

Arrays of exact numbers (new faster methods in Mathematica 8):

LinearSolve[(({{0, 2, 9, 0}, {1, 8, 7, 0}, {8, 2, 0, 0}, {5, 0, 0, 1}})),(({{5}, {3}, {3}, {3}}))]

Machine precision arrays:

LinearSolve[RandomReal[1,{20,20}],RandomReal[1,{20}]]

Linear Algebra: Operations

All standard linear algebra operations, such as:

LeastSquares[RandomReal[1,{100,20}],RandomReal[1,{100}]]

Eigenvalues[RandomReal[1,{10,10}]]

Linear Algebra: Sparse Arrays

Fast, efficient representation of sparse linear equations.

SparseArray[{{4,2}|{3,6}->4,{i_,i_}->2},{8,8}]

MatrixPlot[%]

Linear Algebra: Sparse Arrays

Example: Discretization of 2D heat equation leads to linear equations:

Linear equations as a sparse array:

n=30;

eqns=N[SparseArray[{{i_,j_,i_,j_}->-4,{i_,j1_,i_,j2_}/;Abs[j1-j2]==1->1,{i1_,j_,i2_,j_}/;Abs[i1-i2]==1->1},{n,n,n,n}]]

eqns=Flatten[eqns,{{1,2},{3,4}}]

Overall structure of the linear equations:

MatrixPlot[eqns]

Heat source:

source=ImageData[ColorConvert[ImageResize[ImageCrop[Rasterize["8",ImageSize->n]],{n,n}],"Grayscale"]]-1;

ArrayPlot[source,ImageSize->Small]

Solve discretized heat equations:

sol=LinearSolve[eqns,Flatten[source]];

ListPlot3D[Partition[sol,n]]

Specifying Methods

For when you want to test or compare methods, or you already know the best method to use.

Explore adaptive Monte-Carlo integration:

NIntegrate[Boole[x^2+y^2<1],{x,-1,1},{y,-1,1},Method->"AdaptiveMonteCarlo",PrecisionGoal->2]

Specifying Methods

Compare iterative and direct methods for a large sparse system of equations:

m=SparseArray[{Band[{1,1}]->RandomReal[{1,1},10^5],Band[{1,2}]->2,Band[{2,1}]->2},{10^5,10^5}]

MatrixPlot[Take[m,30,30]]

b=RandomReal[1,10^5];

Iterative method:

Timing[LinearSolve[m,b,Method->"Krylov"];]

Direct method:

Timing[LinearSolve[m,b,Method->"Multifrontal"];]

Specifying Methods

Use a method that switches between algorithms for stiff and non-stiff systems:

Clear [x,y]

{x[t],y[t]}/.First[NDSolve[{x'[t]==y[t],((3 y'[t])/1000)==-x[t]+(1-x[t]^2)y[t]+Sin[10t],x[0]==2,y[0]==0},{x,y},{t,0,10},Method->"StiffnessSwitching"]]

Plot[%,{t,0,10}]

Advanced Method Options: Submethods and Options

Specify submethods and options. See Advanced Documentation.

Submethods: Overall subdivision strategy and specific integration rule:

NIntegrate[Sin[10Sin[x]x],{x,0,10},Method->{"LocalAdaptive",Method->"ClenshawCurtisRule"}]

Options: Search for a minimum starting from a large number of points:

NMinimize[(E^(Sin[50 x]))+(1/4) ((x^2)+(y^2))+Sin[60 (E^y)]-Sin[10 (x+y)]+Sin[70 Sin[x]]+Sin[Sin[80 y]],{x,y},Method->{"RandomSearch","SearchPoints"->250}]

Advanced Method Options: Special Methods

Some methods can be used to do additional work or get results in a different form.

Example: Landing point of a projectile launched over terrain with height h(x):

h[x_]:=5+Cos[x/20]+((3Cos[x/3])/(2+Sin[x/10]));

terrain=Plot[h[x],{x,0,80},Filling->0,PlotRange->{0,30},AspectRatio->Automatic,ImageSize->500]

Projectile under drag, launched with speed s at angle θ:

eqn[s_,θ_]:={{p1[0],p2[0]}=={0,h[0]},{p1'[0],p2'[0]}==s{Cos[θ],Sin[θ]},{p1''[t],p2''[t]}==force[{p1'[t],p2'[t]}]}

force[v_List]:={0,-9.8}-0.02Norm[v]v

The special “EventLocator” method of NDSolve can detect the collision with the ground:

{p1[t],p2[t]}/.First[NDSolve[eqn[30,45°],{p1[t],p2[t]},{t,0,10},Method->{"EventLocator","Event":>p2[t]==h[p1[t]],"EventAction":>Throw[end=t,"StopIntegration"]}]]

Show[terrain,ParametricPlot[%,{t,0,end},PlotStyle->Orange]]

Do the same thing inside a Manipulate:

Manipulate[Show[terrain,ParametricPlot[Evaluate[{p1[t],p2[t]}/.First[NDSolve[eqn[s,θ °],{p1[t],p2[t]},{t,0,10},Method->{"EventLocator","Event":>p2[t]==h[p1[t]],"EventAction":>Throw[end=t,"StopIntegration"]}]]],{t,0,end},PlotStyle->Orange]],{{s,40},1,60,Appearance->"Labeled"},{{θ,45},0,90,Appearance->"Labeled"}]

How to Watch Algorithms: EvaluationMonitor

EvaluationMonitor is a quick and easy way to find out what an algorithm is doing:

FindRoot[Cos[x]==x,{x,0},EvaluationMonitor:>Print[x]]

How to Watch Algorithms: EvaluationMonitor

Watch NIntegrate converge on the more rapidly varying region of this integrand:

Plot[(1/(1+Exp[-x])),{x,-10,20}]

Sow and Reap to get a list of evaluation points:

ListPlot[Last@Reap[NIntegrate[(1/(1+Exp[-3x])),{x,-10,20},EvaluationMonitor:>Sow[x]]]]

Another way to visualize the same points:

ListPlot[Last@Reap[NIntegrate[(1/(1+Exp[-3x])),{x,-10,20},EvaluationMonitor:>Sow[{x,(1/(1+Exp[-x]))}]]],Filling->Axis]

How to Watch Algorithms: StepMonitor

StepMonitor lets you monitor the progress after each algorithm step.

FindMinimum[(1-x)^2+100(-x^2-y)^2+1,{{x,-1.2},{y,1}},StepMonitor:>Print[{x,y}]]

How to Watch Algorithms: StepMonitor

You can access the solution computed so far inside StepMonitor:

Dynamic[Plot[sol,{x,-10,10},PlotRange->{-5,15},Filling->Axis]]

NDSolve[{(∂_(t,t))u[t,x]==(∂_(x,x))u[t,x]+Sin[u[t,x]] ,u[0,x]==(E^(-(x^2))),Derivative[1,0][u][0,x]==0,u[t,-10]==u[t,10]},u,{t,0,50},{x,-10,10},StepMonitor:>(sol=u[t,x])]

Combine with Monitor for a progress indicator:

Monitor[NDSolve[{(∂_(t,t))u[t,x]==(∂_(x,x))u[t,x]+Sin[u[t,x]] ,u[0,x]==(E^(-(x^2))),Derivative[1,0][u][0,x]==0,u[t,-10]==u[t,10]},u,{t,0,50},{x,-10,10},StepMonitor:>(sol=u[t,x];progress=t)],ProgressIndicator[progress,{0,50}]]

Exact and Approximate Numbers: Three Levels of Precision

Exact and Approximate Numbers: Three Levels of Precision

Built-in operations respect the precision of their inputs.

Exact input ⟶ exact output

Sqrt[8]

Precision[%]

Arbitrary precision input ⟶ arbitrary precision output

Sqrt[8`50]

Precision[%]

Machine precision input ⟶ machine precision output

Sqrt[8.0]

Precision[%]

Exact and Approximate Numbers: Combining Precisions

When different precisions are combined, a justifiable output precision is determined:

8`50+2`25

Precision[%]

If any part is MachinePrecision, the result will be MachinePrecision too:

Sin[5`50]+(Sqrt[6.5])

Precision[%]

Exact and Approximate Numbers: Precision & Accuracy

Precision is the number of significant figures. Accuracy is the number of decimal places.

10 digits of precision, low accuracy:

x=10000000.000`10;

{Precision[x],Accuracy[x]}

10 digits of precision, high accuracy:

x=0.0000001`10;

{Precision[x],Accuracy[x]}

Precision & Accuracy Goals in Numerics Functions

Functions try to get enough correct digits (PrecisionGoal), or enough correct digits after the decimal point (AccuracyGoal).

Get either 7 correct significant figures or 2 correct decimal places:

NIntegrate[Exp[10Sin[x]],{x,0,1},PrecisionGoal->7,AccuracyGoal->2]

Precision & Accuracy Goals: Relative Error

When you require a certain number of significant figures, use only PrecisionGoal.

FindRoot[(Sqrt[x])==2/x,{x,1,2π},PrecisionGoal->5,AccuracyGoal->∞]

Precision & Accuracy Goals in Numerics Functions: Zeroes

When the correct answer is zero, a PrecisionGoal can never be satisfied. Include an AccuracyGoal.

NIntegrate[Sin[x],{x,0,2π},AccuracyGoal->10,PrecisionGoal->10]

Using Numerics Functions Together

Often you want to use the result of a numerical function as the input to another function.

Numerically solve a differential equation that depends on a parameter a:

Clear[sol,a,x,y]

sol[a_]/;NumericQ[a]:=y/.First[NDSolve[{y''[x]==-y[x](1+Cos[x+y[x]]),y[0]==1,y'[0]==a},y,{x,0,10}]]

sol[5]

Plot it for different values of a:

Animate[Plot[Evaluate[sol[a][x]],{x,0,10},PlotRange->4],{a,0,1}]

Compute the arc length of the solution as a function of a, using numerical integration:

len[a_]/;NumericQ[a]:=NIntegrate[(Sqrt[1+sol[a]'[x]^2]),{x,0,10}]

Roughly plot the arc length as a function of a:

Plot[len[a],{a,0,1},MaxRecursion->0,PlotPoints->10]

Find the solution with minimum arc length near a==0.8:

min=FindMinimum[len[a],{a,0.8}]

Using Numerics Functions Together: Interpolation

Create fast, pre-computed functions.

Density ρ(r) of a spherically symmetric object recovered from surface density s(r):

s[r_]:=10^-(r^0.3)

ρ[r_]/;r>0:=-(1/π)NIntegrate[((D[s[t],t])/(Sqrt[(t^2)-(r^2)])),{t,r,∞}]

ρ(r) is relatively expensive function:

Timing[Plot[ρ[r],{r,1,100}]]

Pre-interpolate ρ(r) over a range of values:

ρfunc=FunctionInterpolation[ρ[r],{r,1,100},MaxRecursion->10]

Timing[Plot[ρfunc[r],{r,1,100}]]

The interpolation is fast enough to use in other functions:

mass[r_]/;r>0:=NIntegrate[4π (t^2)ρfunc[t],{t,0,r}]

Plot[mass[r],{r,0,100}]

Fast Built-In Functions

Often the shortest way to code something is also the fastest.

data=N[Range[(10^6)]];

Manual loop—slow:

Timing[For[i=1,i≤(10^6),i++,Sin[data[[i]]]];]

Built-in functional operation—fast:

Timing[Map[Sin, data];]

Listable built-in operation—fastest:

Timing[Sin[data];]

Numbers and Arrays: Use Approximate Numbers Early

Exact inputs give exact results, approximate inputs give approximate results:

Clear[f]

f[x_]:=Sin[x]+x

Nest[f,5,2]

Get the corresponding machine precision number:

N[Nest[f,5,2]]

Nest[f,5.0,2]

You may need the exact expression:

N[Nest[f,5,2],50]

D[Nest[f,x,2],x]

If you only need the final approximate number, use approximate numbers from the start:

N[Nest[f,5,20]]//Timing

Nest[f,5.,20]//Timing

Numbers and Arrays: Packed Arrays

Lists and arrays of machine precision numbers (real, complex, integer) may be stored internally as “packed arrays”.

Packed arrays are created automatically by Mathematica:

packed=RandomReal[1,{10^7}];

unpacked=Developer`FromPackedArray[packed];

Packed arrays use less memory:

ByteCount[packed]

ByteCount[unpacked]

Packed arrays are faster:

Sin[packed];//Timing

Sin[unpacked];//Timing

Numbers and Arrays: Packed Arrays

To optimize an algorithm that works with lists of machine numbers, keep everything packed.

v=RandomReal[1,{10^6}];

Issue extra messages about packed arrays:

On["Packing"];

A quantity such as an exact 0 causes this machine precision array to be unpacked:

v[[3]]=0;

Exp[v];//Timing

With a machine precision zero Mathematica can repack the array:

v[[3]]=0.;

Exp[v];//Timing

Switch extra messages back off:

Off["Packing"];

Compilation

Create a fast version of a function that only works on machine numbers.

f=
Compile[{{c,_Complex}},Module[{i=1},FixedPoint[(i++;#^2+c)&,0,50,SameTest->(Re[#]^2+Im[#]^2≥4&)];i],
	CompilationTarget->"C",
	RuntimeAttributes->{Listable}];

Compile produced virtual machine instructions and a native code dynamic library:

InputForm[f]

The compiled function is fast enough to visualize immediately:

ArrayPlot[f[Table[x+I y,{x,-2.1,0.7,0.005},{y,-1.3,1.3,0.005}]],ColorFunction->"GreenPinkTones"]

Zoom in (interactive version: see the CompilationTarget examples):

ArrayPlot[f[Table[x+I y,{x,0,0.5,0.001},{y,0.3,1,0.001}]],ColorFunction->"GreenPinkTones"]

Compilation: Automation

You rarely need to compile because Mathematica does it automatically in many places.

Table[x^2+Cos[x^2],{x,0,1,0.001}];

Plot[(1/(2+Sin[x])),{x,0,5}]

NDSolve[{Sin[x]+y'[x]==Cos[y[x]],y[0]==0},y,{x,0,3}];

Mathematica can’t detect and compile a “hidden” definition like this:

Clear[f];

f[x_]/;NumberQ[x]=Integrate[(1/(1+x^4)),x]

Table[f[x],{x,0,1,0.0001}];//Timing

It can be faster to compile it yourself:

compiled=Compile[{x},Evaluate[Integrate[(1/(1+x^4)),x]]];

g[x_]/;NumberQ[x]:=compiled[x]

Table[g[x],{x,0,1,0.0001}];//Timing

Initialization

m8red[a_]:=RGBColor[0.811765,0.117647,0.145098,a];

Needs["CCodeGenerator`"]

PicturePrint[expr_,opts:OptionsPattern[Cell]]:=CellPrint[TextCell[expr,"Picture",opts,ShowStringCharacters->False]]

OverviewMouseover[header_,subs_]:=CellPrint@Cell[BoxData[ToBoxes@Mouseover[Framed[Style[header,"OverviewSection"],FrameMargins->{{15,15},{5,5}},BoxFrame->2,FrameStyle->None,RoundingRadius->5],

Evaluate@If[subs=="",Framed[Style[header,"OverviewSection",White],FrameMargins->{{15,15},{5,5}},BoxFrame->2,FrameStyle->None,RoundingRadius->5,Background->m8red[1]],Row[{Framed[Style[header,"OverviewSection",White],FrameMargins->{{15,15},{5,5}},FrameStyle->None,RoundingRadius->5,Background->m8red[1]]," » ",Framed[Style[subs,"OverviewSection"],FrameMargins->{{15,15},{5,5}},FrameStyle->m8red[1],BoxFrame->2,RoundingRadius->5,Background->White]}]]]],"OverviewSection"]

LargeOverviewMouseover[header_,subs_]:=CellPrint@Cell[BoxData[ToBoxes@Mouseover[Framed[Style[header,"LargeOverviewSection"],FrameMargins->{{15,15},{5,5}},BoxFrame->2,FrameStyle->None,RoundingRadius->5],

Evaluate@If[subs=="",Framed[Style[header,"LargeOverviewSection",White],FrameMargins->{{15,15},{5,5}},BoxFrame->2,FrameStyle->None,RoundingRadius->5,Background->m8red[1]],Row[{Framed[Style[header,"LargeOverviewSection",White],FrameMargins->{{15,15},{5,5}},FrameStyle->None,RoundingRadius->5,Background->m8red[1]]," » ",Framed[Style[subs,"LargeOverviewSection"],FrameMargins->{{15,15},{5,5}},FrameStyle->m8red[1],BoxFrame->2,RoundingRadius->5,Background->White]}]]]],"LargeOverviewSection"]

MapThread[LargeOverviewMouseover,Transpose[{{"Numerical Functions","Overview and applications"},{"Numerical Methods","Detailed algorithm control"},{"Numbers","Precision and accuracy"},{"Techniques","Programming and performance"}}]];

CreateSpecialCell[expr_]:=CellPrint@Cell[BoxData@ToBoxes@expr,"Picture",CellAutoOverwrite->False]

CreateSpecialTraditionalFormCell[expr_]:=CellPrint@Cell[BoxData@ToBoxes[expr,TraditionalForm],"Picture",CellAutoOverwrite->False]

Block[{y={{0,1.01},{0.2,1.03},{0.4,0.93},{0.6,0.85},{0.8,0.72},{1,0.57}}},ListLinePlot[y,Filling->Axis,InterpolationOrder->3,Mesh->Full,PlotStyle->PointSize[Large]]//CreateSpecialCell]

Block[{path,f,z},
path[t_]:=2+Exp[I t];
f[z_]:=((Sin[z])/((z-2)^3));
Show[Plot3D[Abs[f[x+I y]],{x,0.5,3.5},{y,-1.5,1.5},ColorFunctionScaling->False,ColorFunction->(Darker@Hue[Rescale[Arg[((Sin[#+I #2])/((#+I #2-2)^3))],{-π,π},{0,1}]]&),PlotRange->{Automatic,5},PlotPoints->50,PlotLabel->Style[HoldForm@((Sin[z])/((z-2)^3)),14],AxesLabel->{Re[z],Im[z],Abs[z]},LabelStyle->14],ParametricPlot3D[Append[Through[{Re,Im}[path[t]]],Abs[f[path[t]]]],{t,0,2π},PlotStyle->{Dashing[{}],Yellow,AbsoluteThickness[5]}]]]//CreateSpecialCell

Plot[{Cos[Cos[x]],x},{x,0,1},Epilog->{Red,PointSize[0.03],Point[{0.7390851332151608`,0.7390851332151608`}]}]//CreateSpecialCell

TraditionalForm[Style[HoldForm[f[x]==Integrate[Sin[Exp[t]],{t,0,x}]=="?"],Large]]//CreateSpecialTraditionalFormCell

Plot[Exp[-x^2] Cos[5Exp[x]],{x,-2,2},Epilog->{Red,PointSize[Large],Point[Reverse@MapAt[Composition[Last,Last],FindMinimum[Exp[-x^2] Cos[5Exp[x]],{x,1}],2]]}]//CreateSpecialCell

Plot[Exp[-x^2] Cos[5Exp[x]],{x,-2,2},Epilog->{Red,PointSize[Large],Point[Reverse@MapAt[Composition[Last,Last],NMinimize[{Exp[-x^2] Cos[5Exp[x]]},x],2]]}]//CreateSpecialCell

Block[{z,xy},
{z,xy}=NMinimize[{Sin[x]Cos[2y+1],x^4+y^2<3},{x,y}];
p=Show[Plot3D[Sin[x]Cos[2y+1],{x,-2,2},{y,-2,2},RegionFunction->({x,y}|->x^4+y^2<3),MeshFunctions->{#3&},AxesLabel->{x,y}],Graphics3D[{Red,PointSize[0.05],Point[{x,y,z}/.xy]}]]
]//CreateSpecialCell

Block[{pts=Last[Last[Reap[NIntegrate[Boole[x^2+y^2<1],{x,-1,1},{y,-1,1},Method->"AdaptiveMonteCarlo",EvaluationMonitor:>Sow[{x,y}],PrecisionGoal->2]]]],in,out},
in=Select[pts,(Norm[#]<1)&];
out=Complement[pts,in];RegionPlot[x^2+y^2<1,{x,-1,1},{y,-1,1},Epilog->{Red,Point[in],Lighter[Gray],Point[out]},PlotLabel->Style[x^2+y^2<1,16]]]//CreateSpecialCell

Clear[q,x,y,u]
CreateSpecialTraditionalFormCell[Style[Row[{HoldForm[D[u,{x,2}]+D[u,{y,2}]==q],"⟹",HoldForm[((Subscript[u, i-1,j])-2(Subscript[u, i,j])+(Subscript[u, i+1,j]))+((Subscript[u, i,j-1])-2(Subscript[u, i,j])+(Subscript[u, i,j+1]))==q]},Spacer[30]],20]]

Graphics[{Opacity[0.3],Orange,Rectangle[{2,2}],{Opacity[1],Black,Inset[Style[TraditionalForm[-4(Subscript[u, i,j])],18],{2.5,2.5}]},Blue,Table[Rectangle[p],{p,{{1,2},{3,2},{2,1},{2,3}}}],{Opacity[1],Black,{Inset[Style[TraditionalForm[HoldForm[+(Subscript[u, i-1,j])]],18],{1.5`,2.5`}],Inset[Style[TraditionalForm[+(Subscript[u, i+1,j])//HoldForm],20],{3.5`,2.5`}],Inset[Style[TraditionalForm[+(Subscript[u, i,j-1])//HoldForm],20],{2.5`,1.5`}],Inset[Style[TraditionalForm[+(Subscript[u, i,j+1])//HoldForm],20],{2.5`,3.5`}]}}},PlotRange->{{0,6},{0,6}},ImageSize->{400,400},GridLines->{Range[0,6],Range[0,6]}]//CreateSpecialCell

Block[{x,y,z,xy,p},
{z,xy}=NMinimize[(E^(Sin[50 x]))+(1/4) ((x^2)+(y^2))+Sin[60 (E^y)]-Sin[10 (x+y)]+Sin[70 Sin[x]]+Sin[Sin[80 y]],{x,y},Method->{"RandomSearch","SearchPoints"->250}];
p=Show[Plot3D[(E^(Sin[50 x]))+(1/4) ((x^2)+(y^2))+Sin[60 (E^y)]-Sin[10 (x+y)]+Sin[70 Sin[x]]+Sin[Sin[80 y]],{x,-0.4,0.4},{y,-0.8,0},PlotPoints->50],Graphics3D[{Red,PointSize[0.05],Point[{x,y,z}/.xy]}]]
]//CreateSpecialCell

Block[{data=Block[{n=1000},RandomVariate[ProbabilityDistribution[(1/(r+1)),{r,0,10}],n]Normalize/@RandomVariate[NormalDistribution[],{n,3}]]},Row[{Graphics3D[{Opacity[0.5],White,Point[data]},ViewPoint->{0,0,10},ViewVertical->{0,1,0},Background->Black,ImageSize->300,SphericalRegion->True],SmoothDensityHistogram[Take[data,All,2],ImageSize->300]},Spacer[20]]]//CreateSpecialCell

CreatePictureGridCell[exprs_List,height_]:=CreateSpecialCell[Style[Grid[Map[Image[ImageResize[#,{Automatic,height}],ImageSize->All]&,exprs,{2}],Spacings->{1,1}]]]

CreatePictureGridCell[{{Block[{path,f,z},
path[t_]:=2+Exp[I t];
f[z_]:=((Sin[z])/((z-2)^3));
Show[Plot3D[Abs[f[x+I y]],{x,0.5,3.5},{y,-1.5,1.5},ColorFunctionScaling->False,ColorFunction->(Darker@Hue[Rescale[Arg[((Sin[#+I #2])/((#+I #2-2)^3))],{-π,π},{0,1}]]&),PlotRange->{Automatic,5},PlotPoints->50,LabelStyle->14],ParametricPlot3D[Append[Through[{Re,Im}[path[t]]],Abs[f[path[t]]]],{t,0,2π},PlotStyle->{Dashing[{}],Yellow,AbsoluteThickness[5]}]]],Block[{z,xy},
{z,xy}=NMinimize[{Sin[x]Cos[2y+1],x^4+y^2<3},{x,y}];
p=Show[Plot3D[Sin[x]Cos[2y+1],{x,-2,2},{y,-2,2},RegionFunction->({x,y}|->x^4+y^2<3),MeshFunctions->{#3&},AxesLabel->{x,y}],Graphics3D[{Red,PointSize[0.05],Point[{x,y,z}/.xy]}]]
],},{Block[{x,y,z},{x,y,z}={(2+Cos[v])Cos[u],(2+Cos[v])Sin[u],Sin[v](5/4-Cos[u])};
ParametricPlot3D[{x,y,z},{u,0,2π},{v,0,2π}]],Plot[Evaluate[{x[t],x'[t]}/.First[NDSolve[{x'[t]==x[t](x[t-Pi]-x'[t-1]),x[t/;t≤0]==Cos[t]},x,{t,0,5}]]],{t,0,5}],MatrixPlot[SparseArray[{{4,2}|{3,6}->4,{i_,i_}->2},{8,8}]]}},150]

CreatePictureGridCell[{{,},{,}},150]

CreatePictureGridCell[{{,}},250]

$Line=0;