felipe zapata | 17 Apr 00:36 2013
Picon

Yampa integral function

Hello,
I'm trying to write some simulator of 
physical systems using Yampa.
The fundamental idea is to integrate 
the position and velocity, using
the following algorithm:

1. get the velocity using the formula,

v = v0 + integral ( 0.5 * force * recip mass)

where v0 is the previous velocity.

2. calculate the new position as follows,

x = x0 + integral (v)

and some other steps that are not relevant here.
Therefore I wrote the following piece of 
code:

type Pos   = Double
type Vel   = Double
type Mass  = Double
type Force = Double

integralCoordSF :: Force -> Mass ->  SF (Pos,Vel) (Pos,Vel)
integralCoordSF  force mass = proc (x0,v0) -> do
                                 v <- integralVelSF force mass -< v0
                                 vdt2 <- integral -< v
                                 let x = x0 + vdt2
                                 returnA -< (x,v)
       
integralVelSF :: Force -> Mass ->  SF Double Double       
integralVelSF force mass = proc v0 -> do
                             t1 <- integral -< 0.5 * (force/mass)
                             returnA -< v0 + t1


Now, let's the speed (velocity), force and the mass equal 1.0.  And
the initial potion equal zero. Then using the "embed" function we get that,

*Main>  embed (integralCoordSF  1.0 1.0) ((0.0,1.0),[(0.1,Nothing),(0.1,Nothing)])
[(0.0,1.0),(0.1,1.05),(0.20500000000000002,1.1)]

According to the algorithm, if  dt = 0.1 then 

v = 1.0 + integral (0.5 * 1.0 * recip 1.0) = 1.05

in agreement with the algorithm.  But let's see about the position

x = 0 + integral (v) 

Using the formula of step 2, the position should be equal to 0.105,
but the result given by embed is 0.1.

As a beginner certainly and missing a lot of concepts and details,
can someone please explain me what's going on here?

Thank you for your help!

Felipe Z

 


_______________________________________________
Haskell-Cafe mailing list
Haskell-Cafe <at> haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe
Jerzy Karczmarczuk | 17 Apr 01:48 2013
Picon

Re: Yampa integral function

Le 17/04/2013 00:36, felipe zapata a écrit :
Hello,
I'm trying to write some simulator of 
physical systems using Yampa.
The fundamental idea is to integrate 
the position and velocity, using
the following algorithm:
... ...
Nothing to do with Haskell.
Where do you get this factor 0.5 in your formula:  v = v0 + integral ( 0.5 * force * recip mass) ?
With constant acceleration v=v0+a*Dt =>  1.01, not 1.05

What is that:  let x = x0 + vdt2 ??
Again, a constant acceleration movement gives : x = x0 + v0*Dt + a*Dt^2/2.

Yampa or not Yampa, I suspect that your implementation of physics is simply wrong.
==

In general, even correcting all, you might have reasonable results in some trivial cases, but in general the extrapolating Euler schema is unstable, produces growing errors (e. g. in the oscillating case).

Jerzy Karczmarczuk





_______________________________________________
Haskell-Cafe mailing list
Haskell-Cafe <at> haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe
Jerzy Karczmarczuk | 17 Apr 01:55 2013
Picon

Re: Yampa integral function

Le 17/04/2013 01:48, Jerzy Karczmarczuk a écrit :
> With constant acceleration v=v0+a*Dt =>  1.01, not 1.05
Gosh, trivial errors seem to be contagious. Of course I meant 1.1, not 
1.01. Sorry. JK
L Corbijn | 17 Apr 10:53 2013
Picon

Re: Yampa integral function

Warning: no knowledge of Yampa.

It seems that the integration for the position is always using the velocity from the previous step. Looking at the documentation in source code of Yampa this seems plausible as there is also a function called imIntegrate with the comment

-- "immediate" integration (using the function's value at the current time)

Which I assume would yield your answer. As I don't know Yampa I can only speculate about the reason for this difference, the two things I can think of now are better speed and to prevent trouble (in cases like x = integral y, y = integral x).

Lars


On Wed, Apr 17, 2013 at 1:55 AM, Jerzy Karczmarczuk <jerzy.karczmarczuk <at> unicaen.fr> wrote:
Le 17/04/2013 01:48, Jerzy Karczmarczuk a écrit :

With constant acceleration v=v0+a*Dt =>  1.01, not 1.05
Gosh, trivial errors seem to be contagious. Of course I meant 1.1, not 1.01. Sorry. JK


_______________________________________________
Haskell-Cafe mailing list
Haskell-Cafe <at> haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

_______________________________________________
Haskell-Cafe mailing list
Haskell-Cafe <at> haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

Gmane