PolyMathOrg / PolyMath

Scientific Computing with Pharo
MIT License
169 stars 41 forks source link

Add more PolyMath demos #44

Open SergeStinckwich opened 6 years ago

SergeStinckwich commented 6 years ago

What kind of nice demos we can do with PM 1.0 ?

At the moment, we have the following ones :

olekscode commented 6 years ago

Hello! We could create some demos during the meetings of PharoClub. This would be a good educational exercise. But could you show some examples of how these demos should look like? Maybe the links to the existing ones?

SergeStinckwich commented 6 years ago

Demo about Pi computation

| a b c k pi oldpi oldExpo expo nBits str |
nBits := 10.

a := PMArbitraryPrecisionFloat one asArbitraryPrecisionFloatNumBits: nBits + 16.
b := (a timesTwoPower: 1) sqrt reciprocal.
c := a timesTwoPower: -1.
k := 1.
oldpi := Float pi.
oldExpo := 2.

[| am gm a2 |
    am := a + b timesTwoPower: -1.
    gm := (a * b) sqrt.
    a := am.
    b := gm.
    a2 := a squared.
    c inPlaceSubtract: (a2 - b squared timesTwoPower: k).
    pi := (a2 timesTwoPower: 1) / c.
    expo := (oldpi - pi) exponent.
    expo isZero or: [expo > oldExpo or: [expo < (-1 - nBits)]]] 
            whileFalse: 
                [oldpi := pi.
                oldExpo := expo.
                k := k + 1].
pi asArbitraryPrecisionFloatNumBits: nBits.
str := (String new: 32) writeStream.
pi absPrintExactlyOn: str base: 10.
str contents
SergeStinckwich commented 6 years ago

Lorentz Attractor

|solver ds stepper system dt sigma r b state values diag|
sigma := 10.0.
r := 28.
b := 8.0/3.0.
dt := 0.01.
system := PMExplicitSystem
    block:[:x :t| |c |
        c := Array new:3.
        c at: 1 put: sigma * ((x at: 2) - (x at: 1)).
        c at: 2 put: r* (x at: 1) - (x at: 2) - ((x at:1) * (x at:3)).
        c at: 3 put: (b negated * (x at: 3) + ((x at:1) * (x at:2))).
        c].
stepper := PMRungeKuttaStepper onSystem: system.
solver := (PMExplicitSolver new) stepper: stepper; system: system; dt: dt.
state := #(10.0 10.0 10.0).

values := (0.0 to:100.0 by: dt) collect: [ :t |
    state := stepper doStep: state time: t stepSize: dt].

b := RTGrapher new.
ds := RTData new.
ds dotShape ellipse size: 1; color: (Color red alpha: 0.3).
ds
    points: values;
    x: [:v | [v at: 1] on: Exception do: [ 0.0]];
    y: [:v | [v at: 2] on: Exception do: [ 0.0]].
b add: ds.
b 
AtharvaKhare commented 5 years ago

How about these:

polynomial := PMPolynomial coefficients: #(1 1 1 1).
grapher := RTGrapher new.
points := (-4 to: 4 by: 0.1).

[ polynomial degree >= 0] whileTrue: [ 
    ds := RTData new.
    ds noDot.
    ds label: (polynomial asString).
    ds points: points.
    ds y: polynomial.
    ds x: #yourself.
    ds connectColor: (Color random).
    grapher add: ds.
    polynomial := polynomial derivative.
     ].

grapher legend addText: 'Legend'.
grapher build.
(RTSVGExporter) exportViewAsSVG: (grapher view) filename: 'Successive Derivatives.svg'.

Successive Derivatives

polynomial := PMPolynomial coefficients: #(1 1 1 1).
grapher := RTGrapher new.
points := (-4 to: 4 by: 0.01).
maxDegree := 8.

[ polynomial degree <= maxDegree] whileTrue: [ 
    ds := RTData new.
    ds noDot.
    ds label: (polynomial asString).
    ds points: points.
    ds y: polynomial.
    ds x: #yourself.
    ds connectColor: (Color random).
    grapher add: ds.
    polynomial := polynomial integral.
     ].

grapher legend addText: 'Legend'.
grapher build.
(RTSVGExporter) exportViewAsSVG: (grapher view) filename: 'Successive Integrals.svg'.

Successive Integrals

SergeStinckwich commented 5 years ago

Maybe the implementation of the Batman equation could be a nice demo: https://math.stackexchange.com/questions/54506/is-this-batman-equation-for-real

SergeStinckwich commented 5 years ago

I would like to prepare a demo on Automatic-Differentiation

SergeStinckwich commented 4 years ago

Rumple Polynoms with Fractions

f := [ :x :y |
    (((1335/4) - (x raisedTo: 2)) * (y raisedTo:6)) + ((x raisedTo: 2)*
    ((11 * (x raisedTo: 2)* (y raisedTo: 2)) - (121/1*(y raisedTo: 4)) - 2)) + ((11/2) * (y raisedTo: 8)) + (x/(2*y)) ].
f value: 77617/1 value: 33096/1.
(-54767/66192) asFloat. 
capsulecorplab commented 3 years ago

Has anyone worked on a quaternion estimator demo?