A look at some new features in Wolfram Mathematica

Daniel Lichtblau

Wolfram Research
Champaign IL USA

danl@wolfram.com

Modified version presented by Adam Strzebonski

Kobe University

July 2014

Overview

Over 700 new functions from many areas. These include

I hope to provide a brief glimpse at some of these areas

Acknowledgements

I’m not smart enough to do all this. I thank the many colleagues at Wolfram Research who offered assistance and suggestions in preparing these examples.

Regions (geometric computation)

Region types

Functionality (symbolically and numerically supported)

Region computations

We’ll show some planar regions but there is also support in higher dimensions.

e1 = (x - 2)^2/7 + (y - 1/2)^2/3 - 1;
v2 = {{Sin[Pi/3], Cos[Pi/3]}, {-Cos[Pi/3], Sin[Pi/3]}}.{x - 2, y + 1};
e2 = ({1/8, 1/2}*v2).v2 - 2;
v3 = {{Sin[2 Pi/5], Cos[2 Pi/5]}, {-Cos[2 Pi/5], Sin[2 Pi/5]}}.{x,y + 1};
e3 = ({1/3, 1/7}*v3).v3 - 3/2;

RegionPlot[Evaluate[Thread[{e1,e2,e3}≤0]],{x,-4,7},{y,-5,5}]

regions=Map[ImplicitRegion[#≤0,{x,y}]&,{e1,e2,e3}];
middle=RegionIntersection@@regions;

RegionPlot[regions,PlotRange->All]

mreg=DiscretizeRegion[middle]

RegionNearest[middle,{6.,-3.}]

Some measurements of the region

approxArea=RegionMeasure[N@middle]

AbsoluteTiming[exactArea=RegionMeasure[middle];]

exactArea

exactArea-approxArea

NSolve[{x,y}∈RegionBoundary[regions[[1]]]&&{x,y}∈RegionBoundary[regions[[2]]],{x,y}]

Show[{RegionPlot[Take[regions,2],PlotRange->All],Graphics[{PointSize[Large],Red,Point[{x,y}/.%]}]}]

approxPerimeter=ArcLength[RegionBoundary[mreg]]

Regions as input to solvers

Find the projection of a rotated 4D hypercube on (ℝ^3).

reg=TransformedRegion[Cuboid[{0,0,0,0}],RotationTransform[{{0,0,0,1},{1/2,1/2,1/2,1/2}}]];
ineqs=Reduce[Exists[{t},{x,y,z,t}∈reg],{x,y,z}];
RegionPlot3D[ImplicitRegion[ineqs,{x,y,z}],PlotPoints->100]

Find integer points in a cone/use vector variables.

R=Cone[{{1,2,3},{4,5,6}},3];
Solve[x∈R,x,Integers]

Graphics3D[{{Green,Opacity[0.5],R},{Red,Point[x/.%]}}]

Maximize a function over a 4D ball.

Maximize[{x y+ z t,{x,y,z,t}∈Ball[{0,0,0,0},1]},{x,y,z,t}]

Find a minimum radius disk containing the given 4 points.

SeedRandom[42];
pts=RandomInteger[{-10,10},{4,2}];
{min,pt}=Minimize[{r,pts∈Disk[{a,b},r]},{r,a,b}];
Graphics[{{Green,Disk[{a,b},r]/.pt},{Red,PointSize[Large],Point[pts]}}]

Integrate a function over a 5D ellipsoid.

(∫_({x,y,z,t,u}∈Ellipsoid[{0,0,0,0,0},{1,2,3,4,5}]))((x^2)+(y^2)+(z^2)+(t^2)+(u^2))

Integration over a region of dimension d is performed wrt. the d-dimensional Lebesgue measure.

(∫_(x∈Circle[]))1

NSolve via sparse homotopy

NSolve now supports a homotopy method

This example is around 10 times faster than before. It came from a paper on electronic orbitals by Akihito Kikuchi.

polys={32 s u v (r^4)-336 s u v (r^3)+992 s u v (r^2)-160 s u v r+40 s u v-324 t (u^2) (r^4)+2572 t (u^2) (r^3)-6448 t (u^2) (r^2)+584 t (u^2) r+19936 t (u^2)+156 t (v^2) (r^4)-1068 t (v^2) (r^3)+2248 t (v^2) (r^2)-80 t (v^2) r-32 t (v^2)+26 t ev (r^4)-278 t ev (r^3)+916 t ev (r^2)-126 t ev r-7972 t ev+126 t (r^4)-882 t (r^3)+1748 t (r^2)+1896 t r-12470 t,156 s (u^2) (r^4)-1068 s (u^2) (r^3)+2248 s (u^2) (r^2)-80 s (u^2) r-32 s (u^2)-52 s (v^2) (r^4)+236 s (v^2) (r^3)-32 s (v^2) (r^2)-104 s (v^2) r+48 s (v^2)-26 s ev (r^4)+278 s ev (r^3)-916 s ev (r^2)+126 s ev r-28 s ev-30 s (r^4)+330 s (r^3)-1148 s (r^2)+760 s r-170 s+32 t u v (r^4)-336 t u v (r^3)+992 t u v (r^2)-160 t u v r+40 t u v,156 (s^2) u (r^4)-1068 (s^2) u (r^3)+2248 (s^2) u (r^2)-80 (s^2) u r-32 (s^2) u+32 t v (r^4)-336 s t v (r^3)+992 s t v (r^2)-160 s t v r+40 s t v-324 (t^2) u (r^4)+2572 (t^2) u (r^3)-6448 (t^2) u (r^2)+584 (t^2) u r+19936 (t^2) u+26 u ew (r^4)-278 u ew (r^3)+916 u ew (r^2)-126 u ew r-7972 u ew+126 u (r^4)-882 u (r^3)+1748 u (r^2)+1896 u r-12470 u,-52 (s^2) v (r^4)+236 (s^2) v (r^3)-32 (s^2) v (r^2)-104 (s^2) v r+48 (s^2) v+32 s t u (r^4)-336 s t u (r^3)+992 s t u (r^2)-160 s t u r+40 s t u+156 (t^2) v (r^4)-1068 (t^2) v (r^3)+2248 (t^2) v (r^2)-80 (t^2) v r-32 (t^2) v-26 v ew (r^4)+278 v ew (r^3)-916 v ew (r^2)+126 v ew r-28 v ew-30 v (r^4)+330 v (r^3)-1148 v (r^2)+760 v r-170 v,-13 (s^2) (r^4)+139 (s^2) (r^3)-458 (s^2) (r^2)+63 (s^2) r-14 (s^2)+13 (t^2) (r^4)-139 (t^2) (r^3)+458 (t^2) (r^2)-63 (t^2) r-3986 (t^2)+1000,13 (u^2) (r^4)-139 (u^2) (r^3)+458 (u^2) (r^2)-63 (u^2) r-3986 (u^2)-13 (v^2) (r^4)+139 (v^2) (r^3)-458 (v^2) (r^2)+63 (v^2) r-14 (v^2)+1000,312 (s^2) (u^2) (r^3)-1602 (s^2) (u^2) (r^2)+2248 (s^2) (u^2) r-40 (s^2) (u^2)-104 (s^2) (v^2) (r^3)+354 (s^2) (v^2) (r^2)-32 (s^2) (v^2) r-52 (s^2) (v^2)-52 (s^2) ev (r^3)+417 (s^2) ev (r^2)-916 (s^2) ev r+63 (s^2) ev-60 (s^2) (r^3)+495 (s^2) (r^2)-1148 (s^2) r+380 (s^2)+128 s t u v (r^3)-1008 s t u v (r^2)+1984 s t u v r-160 s t u v-648 (t^2) (u^2) (r^3)+3858 (t^2) (u^2) (r^2)-6448 (t^2) (u^2) r+292 (t^2) (u^2)+312 (t^2) (v^2) (r^3)-1602 (t^2) (v^2) (r^2)+2248 (t^2) (v^2) r-40 (t^2) (v^2)+52 (t^2) ev (r^3)-417 (t^2) ev (r^2)+916 (t^2) ev r-63 (t^2) ev+252 (t^2) (r^3)-1323 (t^2) (r^2)+1748 (t^2) r+948 (t^2)+52 (u^2) ew (r^3)-417 (u^2) ew (r^2)+916 (u^2) ew r-63 (u^2) ew+252 (u^2) (r^3)-1323 (u^2) (r^2)+1748 (u^2) r+948 (u^2)-52 (v^2) ew (r^3)+417 (v^2) ew (r^2)-916 (v^2) ew r+63 (v^2) ew-60 (v^2) (r^3)+495 (v^2) (r^2)-1148 (v^2) r+380 (v^2)+740 (r^3)-3903 (r^2)+7288 r-5102};
vars=Variables[polys];

We are interested in the real-valued solutions.

AbsoluteTiming[hsol = NSolve[polys, vars, Reals];]

NumberForm[Grid[Prepend[Chop[vars/.hsol,10^(-6)],vars],Dividers->All],{6,4}]

Solving a PDE on a region

New numeric PDE  solving capabilities

A simple planar PDE  example

region = 
  RegionDifference[RegionUnion[Disk[], Disk[{5/4,1}]], 
   RegionUnion[Disk[{1/4, 1/4}, 1/5],Rectangle[{1,1},{5/4,5/4}]]];
DiscretizeRegion[region]

Solve a wave equation with initial and boundary conditions

uifWave = Quiet[NDSolveValue[{(∂_(t,t))u[t,x,y] - 
              Inactive[Laplacian][u[t, x, y], {x, y}] == 0, 
          u[0, x, y] == E^(- ((x - 3/2)^2 + (y - 3/2)^2))+3Sin[x]+Cos[y], 
Derivative[1,0,0][u][0, x, y] == 0, 
          DirichletCondition[u[t, x, y] == x+5Sin[y],True]}, 
        u, {t, 0, 2 π}, {x, y} ∈ region]]

Animate[Plot3D[
  uifWave[t, x, y], {x, y} ∈ uifWave["ElementMesh"], 
  PlotRange -> {-5, 10}, Boxed -> False, Axes -> True,
  Mesh -> None], {t, 0, 2 π, 2 π/50}]

DSolve and discrete events

In Mathematica 10 the DSolve function handles discrete events

The classic bouncing ball example

sol = DSolve[{y''[t] == -98/10, y[0] == 135/10, y'[0] == 5, 
   a[0] == 135/10, WhenEvent[y[t] == a[t], y'[t] -> -3/4 y'[t]], 
   WhenEvent[Mod[t, 1] == 0, a[t] -> a[t] - 1]}, {y, a}, {t, 0, 8}, 
  DiscreteVariables -> a]

Plots of the result

Plot[Evaluate[{y[t], a[t]} /. sol], {t, 0, 6}, Filling -> {2 -> 0}, 
 Exclusions -> None, ImageSize -> Medium]

Plot[Evaluate[{y[t], a[t]} /. sol], {t, 0, 6}, Filling -> {2 -> 0}, 
  Exclusions -> None, ImageSize -> Medium,PlotTheme->"Monochrome"]

Plot[Evaluate[{y[t], a[t]} /. sol], {t, 0, 6}, Filling ->{ 2 ->0}, ColorFunction->Function[{x,y},Hue[y]],
  Exclusions -> None, ImageSize -> Medium,FillingStyle->Automatic,PlotTheme->{"NeonColors"}]

Graphs and networks

Many new features

Let’s visit all the capital cities

countries=EntityList["Country"];
caps=EntityValue[EntityList["Country"],"CapitalCity"];
caps=DeleteCases[caps,aa_Missing];

locations=Map[GeoPosition,caps];

Timing[tour=FindShortestTour[locations]]

caps[[tour[[2]]]]

GeoGraphics[GeoPath[caps[[tour[[2]]]]],ImageSize->1000]

Random graph examples

gr=RandomGraph[{20,100}]

cut=FindVertexCut[gr]

Subgraph[gr,cut]

HighlightGraph[gr,Subgraph[gr,Complement[VertexList[gr],cut]],GraphHighlightStyle->"Thick"]

tr=FindSpanningTree[gr]

BreadthFirstScan[tr,1,{"PrevisitVertex"->(Print["Visiting ",#]&)}];

grbig=RandomGraph[{500,4000}]
Timing[cycles5=FindCycle[grbig,{5},All];]
Length[cycles5]

Polynomial algebra

Matrix Smith decomposition (did not make the cut into System context)

??*`PolynomialSmithDecomposition

Example

randomPoly[deg_,max_,x_]:=RandomInteger[{-max,max},deg+1].x^Range[0,deg]
randomMatrix[m_,n_,deg_,max_,x_]:=Table[randomPoly[deg,max,x],{m},{n}]

SeedRandom[1111];
mat=randomMatrix[6,6,8,10,x]

Timing[{uu,ss,vv}=Control`PCS`PolynomialSmithDecomposition[mat,x];]

Check that multiplier matrices are unimodular (free of x).

Expand[{Det[uu],Det[vv]}]//N

Check that middle matrix is diagonal etc.

ss//N

And that we satisfy the appropriate identity.

Expand[uu.mat.vv-ss]

Bigger numerical example

biggermat=randomMatrix[8,10,10,100,x];

Timing[{ubig,sbig,vbig}=Control`PCS`PolynomialSmithDecomposition[N[biggermat,400],x];]

Repeat some of the checks

Chop[Expand[{Det[ubig], Det[vbig]}],10^(-50)]

Chop[Expand[ubig.biggermat.vbig-sbig],10^(-50)]

Gröbner basis conversion matrix (NQRPT)

polys={x^2+3*x*y+y^3-7,x^3*y-5x*y^2-3*x+2*y-4};
{gb,mat}=GroebnerBasis`BasisAndConversionMatrix[polys,{x,y}]

Check result

Expand[mat.polys-gb]

Random processes

Many features

An insurance modeling example from our web pages

Simulate the surplus process for insurance

Insurer's initial surplus = 65

Total annual premium = 61

Number of claims follows a Poisson process with mean = 60

Losses are distributed exponentially with mean = 1

u = 65; c = 61; λ = 60; m = 1;
claims := 
  CompoundPoissonProcess[λ, ExponentialDistribution[1/m]];
InsuranceProcess = 
  TransformedProcess[u + c t - x[t], x \[Distributed] claims, t];

We’ll look at 50 random samples

samples = RandomFunction[InsuranceProcess, {0, 120, 1}, 50]
ListLinePlot[samples, ImageSize -> Medium]

Some of the basic stats on these samples:

Mean[samples[120]]
Min[samples["SliceData",120]]
StandardDeviation[samples[120]]

Inactive and formatting

Mathematica 10 has a function that renders operators inactive (very useful for formatting)

(plus = Table[ 
    Inactive[Plus][a, b] == a + b, {a, 7}, {b, 7}]) // TableForm

Grid[Join[{{Item[Style["Addition Table", FontSize -> 18], 
      Alignment -> Left, Background -> GrayLevel[.9]], SpanFromLeft, 
     SpanFromLeft, SpanFromLeft, SpanFromLeft}}, plus]
  ]

And we can “activate” to evaluate.

Activate[plus]

This is similar to Mathematica’s Hold and Release, but is more “math-aware”. Example: Derivative can be used on an inactivated LaplaceTransform in the way one might expect.

Inactive[LaplaceTransform][a t^2, t, s]

D[%, s]