# Mathematica

## Project !

On Tuesday, towards the end of the exam, we had a bit of an overview something about the project. The project contains also a creative aspect Mathematica. Try to make that part printable. It can be whatever you like. The only requirement is that it should be creative and that it should consist of several graphics objects put together. A simple surface for example would not quality. Also, something like a Snowman or some object which you have found code for on the internet would give no credit.
Here is the Mathematica notebook for 2020.
Some pictures from the mathematica project:

## Mandelbrot set area

```M=Compile[{x,y},Module[{z=x+I y,k=0},
While[Less[Abs[z],2]&&Less[k,1000],z=N[z^2+x+I y];++k];Floor[k/1000]]];
9*Sum[M[-2+3 Random[],-1.5+3 Random[]],{1000000}]/1000000
```

## Checking partial differential equations

Here is an example to verify that a function satisfies the heat equation
 ```f[t_,x_]:=(1/Sqrt[t])*Exp[-x^2/(4t)]; Simplify[D[f[t,x],t] == D[f[t,x],{x,2}]] ```

We did not need FullSimplify here, but here is the analogue stronger nudge. You can also append it to the end after a // this is handy when you want to avoid brackets.
 ```D[f[t,x],t] == D[f[t,x],{x,2}] // FullSimplify ```

When using Mathematica (like an RPN calculator), one can avoid brackets
 8*1231434 // Sqrt // Log // N

## Cobb Douglas by data fitting

In the summary part of the third week, (Problem 300.5) you are asked to visualize the original Cobb-Douglas data points. Here is the code.
 ```(* Historical capital, labor and production data used by Cobb-Douglas *) A = {{1899, 100, 100, 100}, {1900, 107, 105, 101}, {1901, 114, 110, 112}, {1902, 122, 118, 122}, {1903, 131, 123, 124}, {1904, 138, 116, 122}, {1905, 149, 125, 143}, {1906, 163, 133, 152}, {1907, 176, 138, 151}, {1908, 185, 121, 126}, {1909, 198, 140, 155}, {1910, 208, 144, 159}, {1911, 216, 145, 153}, {1912, 226, 152, 177}, {1913, 236, 154, 184}, {1914, 244, 149, 169}, {1915, 266, 154, 189}, {1916, 298, 182, 225}, {1917, 335, 196, 227}, {1918, 366, 200, 223}, {1919, 387, 193, 218}, {1920, 407, 193, 231}, {1921, 417, 147, 179}, {1922, 431, 161, 240}}; B = Table[{A[[k, 3]], A[[k, 2]], A[[k, 4]]}, {k, Length[A]}]; f[x_, y_] := x^(3/4) y^(1/4); (* Cobb- Douglas obtained this by fitting *) S1 = Graphics3D[Table[{Red,Sphere[B[[k]], 10]}, {k, Length[B]}]]; S2 = Plot3D[x^(3/4) y^(1/4),{x, 100, 193}, {y, 100, 400}]; S = Show[{S1, S2},Axes->True, Ticks->None,AxesLabel-> {"C", "L", "P"}] ```

## Curves

Plotting a function
```ParametricPlot[ {Cos[5t],Sin[7t]},{t,0,2Pi}]
ParametricPlot3D[{Cos[5t],Sin[7t],t+Sin[t]},{t,0,2Pi}]
ListPlot[Table[Prime[k],{k,1000}],Joined->True]
```
Velocity and Acceleration of the primes:
```f=Table[Prime[k],{k,10000}];
d[f_]:=Table[f[[k+1]]-f[[k]],{k,Length[f]-1}];
a[f_]:=Table[f[[k+2]]-2*f[[k+1]]+f[[k]],{k,Length[f]-2}];
SmoothHistogram[d[f],PlotRange->All,Filling->Bottom,FillingStyle->Yellow]
SmoothHistogram[a[f],PlotRange->All,Filling->Bottom,FillingStyle->Red]
```

## Surfaces

Contourplots of functions of three variables
```ContourPlot3D[ x^2+y^2+z^2 == 3, {x,-2,2},{y,-2,2},{z,-2,2}]
```
Here were templates for the rare elements surfaces:
```      r=2+Cos[7u];R=6+r*Cos[v];ParametricPlot3D[{R*Cos[u],R*Sin[u],Sin[v]},{u,0,2Pi},{v,0,2Pi}]
r=2+Cos[5u];R=1+r*Cos[v];ParametricPlot3D[{R*Cos[u]*Sin[v],R*Sin[u]*Sin[v],Cos[v]},{u,0,2Pi},{v,0,Pi}]
r=(2+v*Cos[u/2]);R={r*Cos[u],r*Sin[u],-v*Sin[u/2]};ParametricPlot3D[R,{u,0,2 Pi},{v,-1,1}]
r=(2 + Cos[3 t/2]);ParametricPlot3D[{r*Cos[t] + Sin[s],r*Sin[t],Sin[3t/2]},{t,0,4 Pi},{s,0,2Pi}]
ContourPlot3D[((x^2+y^2)^2-x^2+y^2)^2+z^2==0.01,{x,-1.4,1.4},{y,-1,1},{z,-1,1}]
ContourPlot3D[x*y^2+y*z^2+z*x^2==1,{x,-4,4},{y,-4,4},{z,-4,4}]
ContourPlot3D[Abs[x]^(1/3)+Abs[y]^(1/3)+Abs[z]^(1/3)==1.2,{x,-1,1},{y,-1,1},{z,-1,1}]
t=GoldenRatio;w=1; a=1.4;
ContourPlot3D[8(x^2-t^4 y^2)(y^2-t^4 z^2)(z^2-t^4 x^2)(x^4+y^4+z^4-2x^2*y^2-2x^2*z^2-2y^2*z^2)+
(3+5t) (x^2+y^2+z^2-w^2)^2 (x^2+y^2+z^2-(2-t) w^2)^2 w^2==0,{x,-a,a},{y,-a,a},{z,-a,a}]
```
Graphs and contour plots of functions of two variables
```Plot3D[ Sin[x y]/x,{x,-4,4},{y,-4,4}]
ContourPlot[ Sin[x y]/x,{x,-4,4},{y,-4,4}]
```
Tweaking parameters:
```ContourPlot3D[x^2*y+y^2*z + z^2*x == 1,{x,-4,4},{y,-4,4},{z,-4,4},
Extrusion->0.3,Boxed->False,Axes->False,Mesh->False,PlotPoints->40]
```

## Lagrange problems

Solve a Lagrange problem with 2 variables
```f=2x^2+4 x y;  g=x^2 y;
Solve[{D[f,x]==L*D[g,x],D[f,y]==L*D[g,y],g==1},{x,y,L}]
```
With 3 variables:
```f=2x^2+4 x*y+z;     g=x^2 y + z;   c=1;
Solve[{D[f,x]==L*D[g,x], D[f,y]==L*D[g,y], D[f,z]==L*D[g,z], g==c},{x,y,z,L}]
```
With 3 variables and two constraints
```f=x y z;  g=x*y+2 y z+2 x z;  h=2x+2y+4z; d=4; c=1
Solve[{D[f,x]==L*D[g,x] + M*D[h,x],
D[f,y]==L*D[g,y] + M*D[h,y],
D[f,z]==L*D[g,z] + M*D[h,z],
g==c, h==d},{x,y,z,L,M}]
```

## Classifying critical points

Here is example code on how to compute the gradient and the discriminant
```f=x^3 y^3- x y^2;
D[f,{x,2}]*D[f,{y,2}]-D[f,{x,1},{y,1}]^2
```

## Classifying critical points

This produces a nice table.
```f=4 x y - x^3 y - x y^3;
ClassifyCriticalPoints[f_,{x_,y_}]:=Module[{X,P,H,g,d,S}, X={x,y};