上面的是鸟叔骑马舞函数
图像如下

ContourPlot3D[((z - 3 ArcTan[x - 3])^2/(3 (1 + 3 E^(-(2 x/3 - 2)^2))) + (3 + x/10)^2*(y + 1/3 (Sqrt[x^2 + 8] + x) + 1/2)^2/(16 (E^(-x^2/3) + 1)) - 1) ((z - 3 ArcTan[x - 3])^2/(3 (1 + 3 E^(-(2 x/3 - 2)^2))) + (3 + x/10)^2*(y - 1/3 (Sqrt[x^2 + 8] + x) - 1/2)^2/(16 (E^(-x^2/3) + 1)) - 1) ((x - (x + y + z)/9 + 7)^2 + (z - (x + y + z)/9 + 2)^2 + (y - (x + y + z)/9 - 3)^2 - 2) ((x - (x - y + z)/9 + 7)^2 + (z - (x - y + z)/9 + 2)^2 + (y + (x - y + z)/9 + 3)^2 - 2) == 4, {x, -10, 10}, {y, -10, 10}, {z, -10, 10}, ContourStyle -> Pink, Mesh -> None]

图像如下

Nordstrand[x_, y_,
z_] := (2 (4/3 x)^2 + 2 y^2 + z^2 - 1)^3 - (4/3 x)^2 z^3/10 -
y^2 z^3;
Kuska[x_, y_, z_] := (2*x^2 + y^2 + z^2 - 1)^3 - (1/10)*x^2*z^3 -
y^2*z^3;
Taubin[x_, y_, z_] := (x^2 + (3/2)^2 y^2 + z^2 - 1)^3 -
x^2 z^3 - (3/2)^2/20 y^2 z^3;
Trott[x_, y_, z_] :=
320 *((x^2 + (3/2)^2 y^2 + z^2 - 1)^3 - (x^2) z^3 - (3/2)^2/
20 y^2 z^3)
Manipulate[
Switch[type,
"Nordstrand",
ContourPlot3D[
Nordstrand[x, y, z] == 0, {x, -2, 2}, {y, -2, 2}, {z, -2, 2},
Boxed -> False, PlotPoints -> points, MaxRecursion -> max,
PlotRange -> All, ViewPoint -> {2, .1, .5}, BoxRatios -> {1, 1, 1},
Axes -> False, Mesh -> mesh,
ContourStyle -> Directive[col, Specularity[White, 10]],
Epilog ->
Inset[Framed[Style[TraditionalForm[Nordstrand[x, y, z] == 0], 13],
Background -> LightYellow], {Center, Bottom}, {Center,
Bottom}], ImageSize -> {400, 350}]
, "Kuska",
ContourPlot3D[
Kuska[x, y, z] == 0, {x, -0.9, 0.9}, {y, -1.2, 1.2}, {z, -1.2,
1.4}, Boxed -> False, PlotPoints -> points, MaxRecursion -> max,
ViewPoint -> {2, .1, .5}, Axes -> False,
ContourStyle -> Directive[col, Specularity[White, 10]],
Mesh -> mesh,
Epilog ->
Inset[Framed[Style[TraditionalForm[Kuska[x, y, z] == 0], 13],
Background -> LightYellow], {Center, Bottom}, {Center, Bottom}],
ImageSize -> {400, 350}]
, "Taubin",
ContourPlot3D[
Taubin[x, y, z] == 0, {x, -2, 2}, {y, -2, 2}, {z, -2, 2},
Boxed -> False, Axes -> False, PlotPoints -> points,
ContourStyle -> Directive[col, Specularity[White, 10]],
PlotRange -> All, ViewPoint -> {1, 1, .2}
, Epilog ->
Inset[Framed[Style[TraditionalForm[Taubin[x, y, z] == 0], 13],
Background -> LightYellow], {Center, Bottom}, {Center, Bottom}],
ImageSize -> {400, 350}, Mesh -> mesh]
, "Trott",
ContourPlot3D[
Trott[x, y, z] == 0, {x, -3/2, 3/2}, {y, -1, 1}, {z, -3/2, 3/2},
Boxed -> False, PlotPoints -> points, MaxRecursion -> max,
ViewPoint -> {0.26, 3, 0.22}, Axes -> False, Mesh -> mesh,
ContourStyle -> Directive[col, Specularity[White, 10]]
, Epilog ->
Inset[Framed[Style[TraditionalForm[Trott[x, y, z] == 0], 13],
Background -> LightYellow], {Center, Bottom}, {Center, Bottom}],
ImageSize -> {400, 350}]
]
, {{type, "Trott", "heart equation"}, {"Kuska", "Nordstrand",
"Taubin", "Trott"}}
, {{col, Red, "colour"}, Red}
, {{points, 5, "number of sample points"}, 2, 12, 1,
Appearance -> "Labeled"}
, {{max, 1, "max recursion"}, {0, 1, 2, 3, 4, 5},
ControlType -> SetterBar}
, {{mesh, True, "mesh"}, {True, False}}
, ContinuousAction -> False, SaveDefinitions -> True,
TrackedSymbols -> Manipulate
]
图像如下

下午继续

上面的那个函数用我自己电脑上的mathematica画出来的样子

这个代码一大段啊
CenterPointPointArc[c_, p1_, p2_] :=
Module[{a1, a2},
If[(1000 Norm[p2 - p1] > Norm[p2 - c]),
a1 = Mod[ArcTan @@ (p1 - c), 2 \[Pi]];
a2 = Mod[ArcTan @@ (p2 - c), 2 \[Pi]];
If[a1 > a2, a1 -= 2 \[Pi]];
Arc[c, Norm[p1 - c], {a1, a2}], Arc[p1, p2]]]
Render[Arc[c_, r_, {a_, b_}]] :=
If[Abs[a - b] >= 2 \[Pi], Circle[c, r],
Circle[c, r, If[a < b, {a, b}, {b, a}]]]
Render[Arc[a_, b_]] := Line[{a, b}]
Clockwise[radiusVector : {rx_, ry_}, directionVector : {dx_, dy_}] :=
dy rx - dx ry < 0
ClockwiseVector[{x_, y_}] := {y, -x}
CounterclockwiseVector[{x_, y_}] := {-y, x}
TerminalTangentVector[Arc[{xc_, yc_}, r_, {a_, b_}]] :=
If, Sin}],
ClockwiseVector[{Cos, Sin}]]
TerminalTangentVector[Arc[a_, b_]] := b - a
ArcPointArc[arc : Arc[{xc_, yc_}, r_, {a_, b_}], p : {xp_, yp_}] :=
PointTangentPointArc[{xc, yc} + r {Cos, Sin},
TerminalTangentVector[arc], p]
ArcPointArc[arc : Arc[pa_, pb_], p_] :=
PointTangentPointArc[pb, pb - pa, p]
PointTangentPointArc[{xb_, yb_}, tv : {tx_, ty_}, p : {xp_, yp_}] :=
Module[{x, y, aa, bb},(*Solve[{Norm[{xp,yp}-{x,y}]==Norm[{xb,yb}-{x,
y}],{x,y}=={xb,yb}+a {-ty,tx}},{x,y,a}]*)
d = 2 (ty xb - ty xp - tx yb + tx yp);
If[d == 0,
Arc[{xb, yb}, {xp, yp}], {x,
y} = {(ty xb^2 - ty xp^2 - 2 tx xb yb - ty yb^2 + 2 tx xb yp +
2 ty yb yp - ty yp^2)/
d, -((-tx xb^2 + 2 tx xb xp - tx xp^2 - 2 ty xb yb +
2 ty xp yb + tx yb^2 - tx yp^2)/d)};
aa = Mod[ArcTan @@ ({xb, yb} - {x, y}), 2 \[Pi]];
bb = Mod[ArcTan @@ ({xp, yp} - {x, y}), 2 \[Pi]];
Arc[{x, y}, Norm[{xp, yp} - {x, y}],
If[Clockwise[{xb, yb} - {x, y}, tv], {aa,
bb - If, 0]}, {aa,
bb + If, 0]}]]]]
ToPointList[Arc[c_, r_, {a_, b_}], n_: 64] :=
Table[N[c + r {Cos[\[Alpha]], Sin[\[Alpha]]}], {\[Alpha], a,
b, (b - a)/(1 + Abs[IntegerPart[n (b - a)/2 \[Pi]]])}]
ToPointList[Arc[a_, b_], n_: 64] := {a, b}
ReverseArc[Arc[c_, r_, {a1_, a2_}]] := Arc[c, r, {a2, a1}]
ReverseArc[Arc[a_, b_]] := Arc
DiagramArcCircle[Arc[c_, r_, {a_, b_}]] :=
If[Abs[r] > 500,
DiagramArcCircle[
Arc[c + r {Cos[a], Sin[a]}, c + r {Cos, Sin}]], {Circle[c,
r], Arrow[{c, c + r {Cos[a], Sin[a]}}],
Arrow[{c, c + r {Cos, Sin}}]}]
DiagramArcCircle[Arc[a_, b_]] :=
Module[{vab = 10 Normalize[(b - a)],
vabperp}, {vabperp = Reverse[vab] {-1, 1};
Line[{a - vab, b + vab}], Line[{a - vabperp, a + vabperp}],
Line[{b - vabperp, b + vabperp}]
}]
DiagramArc[Arc[c_, r_, {a_, b_}]] :=
If[Abs[r] > 500,
DiagramArc[
Arc[c + r {Cos[a], Sin[a]}, c + r {Cos, Sin}]], {Render[
Arc[c, r, {a, b}]], Arrow[{c, c + r {Cos[a], Sin[a]}}],
Arrow[{c, c + r {Cos, Sin}}]}]
DiagramArc[Arc[a_, b_]] :=
Module[{vab = (b - a), vabperp}, {vabperp = 10 Reverse[vab] {-1, 1};
Line[{a - vab, b + vab}], Line[{a - vabperp, a + vabperp}],
I

Line[{b - vabperp, b + vabperp}]}]
Manipulate[
DynamicModule[{dPts, r, \[Alpha]},
LocatorPane[Dynamic[pts, (dPts = # - pts;
r = Norm[#[[6]] - pts[[3]]];
\[Alpha] =
Mod[(ArcTan @@ (#[[2]] - pts[[3]]) +
ArcTan @@ (#[[4]] - pts[[3]]))/2., 2 \[Pi]];
pts = {{0, #[[1, 2]]}, #[[3]] +
r Normalize[#[[2]] - pts[[3]]], #[[3]], #[[3]] +
r Normalize[#[[4]] - pts[[3]]], {0, #[[5, 2]]}, #[[3]] +
r {Cos[\[Alpha]], Sin[\[Alpha]]}}) &],
Dynamic[Module[{center, bottomArc, middleArc, topArc},
middleArc = CenterPointPointArc[pts[[3]], pts[[2]], pts[[4]]];
bottomArc = ArcPointArc[ReverseArc[middleArc], pts[[1]]];
topArc = ArcPointArc[middleArc, pts[[5]]];
Graphics[{With[{rightHalf =
Join[Reverse@ToPointList, ToPointList[middleArc],
ToPointList[topArc]]}, {Hue[h, .75],
Polygon[Join[Reverse[({-1, 1}*#) & /@ rightHalf],
rightHalf]]}],
If[! scl, {}, {AbsoluteThickness[1], Dashed,
Arrowheads[Medium], Line[{{0, -3.5}, {0, 2.5}}], Black,
DiagramArcCircle[middleArc], Gray, DiagramArc[topArc],
DiagramArc
}
]
}, PlotRange -> {2.5 {-1, 1}, {-3.5, 2.5}}]]],
Appearance ->
If[! sl, None,
Graphics[{Opacity[.35], Disk[]}, ImageSize -> 12,
Background -> None]]]], {{pts, {{0, -2.5}, {2., -0.8}, {1.3,
0.1}, {1.3, 1.3}, {0, 0.8}, {2.4, 0.5}}}, {-2.5`, -3.5`}, {2.5,
2.5}, ControlType -> None}, {{h, .9, "color"}, .8,
1}, {{scl, True, "show construction lines"}, {True, False}}, {{sl,
True, "show control points"}, {True, False}},
SaveDefinitions -> True,
Bookmarks -> {"frame 1" :> {h = 0.9`,
pts = {{0, -2.5`}, {2.`, -0.8`}, {1.3`, 0.1`}, {1.3`, 1.3`}, {0,
0.8`}, {2.4`, 0.5`}}, scl = True, sl = True},
"frame 2" :> {h = 0.9`,
pts = {{0, -2.5`}, {2.023599064419524`,
0.11608691717489805`}, {1.3050000000000002`,
1.04`}, {2.1972423440391537`, 1.797564254372865`}, {0,
0.8`}, {2.469278325693163`, 0.9197669749148628`}}, scl = True,
sl = True},
"frame 3" :> {h = 0.9`,
pts = {{0, -2.5`}, {1.671508330658035`, -0.07160768419360086`}, \
{1.3050000000000002`, 1.04`}, {2.1972423440391546`,
1.7975642543728647`}, {0, 0.8`}, {2.431749417236919`,
0.7230840004728973`}}, scl = True, sl = True},
"frame 4" :> {h = 0.9`,
pts = {{0, -2.5`}, {1.439645555533312`,
0.6316244689319328`}, {1.3050000000000002`,
1.04`}, {1.6327864540426722`, 1.318309253432457`}, {0,
0.8`}, {1.718938206966032`, 0.9235733672489551`}}, scl = True,
sl = True}}]

Plot[1/4 ((3 Sqrt[5] - 5)/8 Abs[x] + (Sqrt[5] + 1)/
2 - (x^2 - ((29 + 8 Sqrt[5])/29)^2)/
Abs[x^2 - ((29 + 8 Sqrt[5])/29)^2] ((5 Sqrt[5] - 3)/8 Abs[x] - (
Sqrt[5] + 1)/2) + (Sqrt[5] - 1) (Abs[x] - 4) +
Sin[5/2 Pi x ((1 + Sqrt[5])/
2 + ((x - 2) (x - 1) (x + 1) (x + 2) (3 - Sqrt[5]))/(
2 Abs[(x - 2) (x - 1) (x + 1) (x + 2)]))] Abs[(
3 Sqrt[5] - 5)/8 Abs[x] + (Sqrt[5] + 1)/
2 - (x^2 - ((29 + 8 Sqrt[5])/29)^2)/
Abs[x^2 - ((29 + 8 Sqrt[5])/29)^2] ((5 Sqrt[5] - 3)/
8 Abs[x] - (Sqrt[5] + 1)/2) - (Sqrt[5] - 1) (Abs[x] -
4)]) (1 - (x^2 - 18)/Abs[x^2 - 16]) +
1/4 ((127 - 95 Sqrt[5])/
29 + (137 Sqrt[5] - 105)/
29 Sin[5 Pi x]) (1 - ((x + 8) (x + 2 (Sqrt[5] + 1)))/
Abs[(x + 8) (x + 2 (Sqrt[5] + 1))]) +
1/4 (-Sqrt[(9 - 3 Sqrt[5])^2 - (x - 11 + Sqrt[5])^2] + 13 -
7 Sqrt[5] +
1/2 (-Sqrt[(3 - Sqrt[5])^2 - (x - 11 + Sqrt[5])^2] + (
388 - 182 Sqrt[5])/
29 - ((x - 8) (x - 14 + 2 Sqrt[5]))/
Abs[(x - 8) (x - 14 +
2 Sqrt[5])] (-Sqrt[(3 - Sqrt[5])^2 - (x - 11 + Sqrt[
5])^2] + (366 - 224 Sqrt[5])/29)) +
Sin[5 Pi x] Abs[-Sqrt[(9 - 3 Sqrt[5])^2 - (x - 11 + Sqrt[
5])^2] + 13 - 7 Sqrt[5] -
1/2 (-Sqrt[(3 - Sqrt[5])^2 - (x - 11 + Sqrt[5])^2] + (
388 - 182 Sqrt[5])/
29 - ((x - 8) (x - 14 + 2 Sqrt[5]))/
Abs[(x - 8) (x - 14 +
2 Sqrt[5])] (-Sqrt[(3 - Sqrt[5])^2 - (x - 11 + Sqrt[
5])^2] + (366 - 224 Sqrt[5])/29))]) (1 - ((x -
2 Sqrt[5] - 2) (x - 20 + 4 Sqrt[5]))/
Abs[(x - 2 Sqrt[5] - 2) (x - 20 + 4 Sqrt[5])]), {x, -9, 12},
PlotRange -> {-10, 11}, AspectRatio -> 1,
ColorFunction -> Function[{x, y}, Hue[y]]]

嗯嗯...这个函数是有用的

还有一个叫生命游戏的东西...
1970年，一位叫Conway的神人，发明了一种叫做“生命”的小游戏：搞一个无限大的正方形网格，每个格子叫做一个细胞，这样每个细胞周围都有八个细胞。细胞有活与死两种状态。每回合每个细胞的死活都按以下的规则变化：
1、若它周围的八个细胞中恰有三个活着，那无论它原来是死是活，这回合它都活着；
2、若它周围的八个细胞中恰有两个活着，那它保持原来的死活状态；
3、其他情况下，无论它原来是死是活，这回合都会死去。
这个游戏可以看做是对生命都某种模拟。周围的细胞少了，会死于寂寞；周围的细胞多了，会死于拥挤。

这就是原图
在mathematica里面用代码
DynamicModule[{a = Binarize[(*把图片复制到这里*)]}, Dynamic[a = MorphologicalTransform[a, "Life"]]]

就会发生这种事情