"Graphical" solve of nonlinear differentialequation with the propeller

Reinhard · 2011-12-13 00:10

Hi,
if someone not interesst in math, skip it ;-))

/*
    Solve the Differential Equation x'(t) = r * x(t) * ( 1 - x(t) )
    with Euler method and send x(t) as Analog out to P0 
    (RC lowpass recommended )
    (P0 observed with an oszi)
    ------------------------------------------------------------------
*/

#include <propeller.h>
#include <math.h>

#define P0 1<<0
#define CRTMODE (1<<28) + (1<<27) + 0 
#define scale  16777216        // 2^32 / 256

typedef float (*Fct)(float x);

void msleep(int t)
{
    waitcnt((CLKFREQ/1000)*t+CNT);
}

float Log_EQ (float x)
{
    float dt = 0.001;
    float r = 0.7;
    return x + dt * r * (1.0 - x);
}

void main ()
{
    int duty;
    Fct eulerstep;

    float x;
    int it;
    
    DIRA |= P0;                    // P0 = OUTPUT
    CTRA = CRTMODE;                // set Counter Mode: DUTY SINGLE ENDED

    while(1)
    {
        x = 0.1;
        eulerstep = Log_EQ;
    
        for (it = 0; it <= 4000; it ++)
        {
            x = eulerstep(x);
            duty = (int)(x * 255);
            FRQA = (duty * scale);        // set FRQA Register 

        }
    FRQA  = 0;
    msleep(100);
    }
}

I think about a alternative solution to do the numerical integration with integer by use the counter ( PSHA += FREQA) property.
Will try it, but now I have go to work .

best regards
Reinhard

ersmith · 2011-12-13 06:48

Very cool demo, Reinhard! Thanks for posting it.

(It's great to see that someone is testing our floating point code, too. :-))

Regards,
Eric

Reinhard · 2011-12-14 12:43

Hi,
here is the fixed point version (much smaller and faster)

/* A Simple C program to illustrate the use of Fixed Point Operations */

/*
    Solve the Differential Equation x'(t) = r * x(t) * ( 1 - x(t) )
    with Euler method and send x(t) as Analog out to P0 
    (RC lowpass recommended )
    (P0 observed with an oszi)

    compile & load with:
    propeller-elf-gcc -mlmm -Os -o euler_fix.elf  euler_fix.c -lm
    propeller-load -pcom7 -r euler_fix.elf  
    ------------------------------------------------------------------
*/

#include <propeller.h>
#include <math.h>

#define P0 1<<0
#define CRTMODE (1<<28) + (1<<27) + 0 


/* DEFINE THE MACROS */
/* The basic operations perfomed on two numbers a and b of fixed
point q format returning the answer in q format */
#define FADD(a,b) ((a)+(b))
#define FSUB(a,b) ((a)-(b))
#define FMUL(a,b,q) (((a)*(b))>>(q))
#define FDIV(a,b,q) (((a)<<(q))/(b))
/* The basic operations where a is of fixed point q format and b is
an integer */
#define FADDI(a,b,q) ((a)+((b)<<(q)))
#define FSUBI(a,b,q) ((a)-((b)<<(q)))
#define FMULI(a,b) ((a)*(b))
#define FDIVI(a,b) ((a)/(b))
/* convert a from q1 format to q2 format */
#define FCONV(a, q1, q2) (((q2)>(q1)) ? (a)<<((q2)-(q1)) : (a)>>((q1)-(q2)))
/* the general operation between a in q1 format and b in q2 format
returning the result in q3 format */
#define FADDG(a,b,q1,q2,q3) (FCONV(a,q1,q3)+FCONV(b,q2,q3))
#define FSUBG(a,b,q1,q2,q3) (FCONV(a,q1,q3)-FCONV(b,q2,q3))
#define FMULG(a,b,q1,q2,q3) FCONV((a)*(b), (q1)+(q2), q3)
#define FDIVG(a,b,q1,q2,q3) (FCONV(a, q1, (q2)+(q3))/(b))
/* convert to and from floating point */
#define TOFIX(d, q) ((int)( (d)*(double)(1<<(q)) ))
#define TOFLT(a, q) ( (double)(a) / (double)(1<<(q)) )



void msleep(int t)
{
    waitcnt((CLKFREQ/1000)*t+CNT);
}


void main ()
{
    int duty;
    float x;
    int it;
    float dt = 0.001;
    float r = 0.7;
    int q = 14;
    int DT,R,X,ONE,T0,T1,T2,T3;
    
    DIRA |= P0;                            // P0 = OUTPUT
    CTRA = CRTMODE;                        // set Counter Mode: DUTY SINGLE ENDED

    while(1)
    {
        x = 0.1;
        DT  = TOFIX(dt,q);
        R   = TOFIX(r,q);
        ONE = TOFIX(1.0,q);
        X   = TOFIX(x,q);
        for (it = 0; it <= 4000; it ++)
        {
            T3 = FSUB(ONE,X);           // perform x = x + dt * r * (1.0 - x)
            T2 = FMUL(R,T3,q);          // with fixed point operations
            T1 = FMUL(DT,T2,q);
            T0 = FADD(X,T1);
            X  = T0;
            duty = X;
            FRQA = (duty << 18);        // set FRQA Register in full range
            msleep(1);                  // sleep a while, avoid be to fast for the RC Lowpass on P0 !
        }
    FRQA  = 0;
    msleep(100);
    }
}

used macros found on ARM website:

Reinhard · 2011-12-14 13:18

\myC\propeller_chip\gcc\euler_schritt>propeller-elf-gcc -mcog -Os -o euler_fix.elf euler_fix.c

\myC\propeller_chip\gcc\euler_schritt>propeller-load -pcom7 -r euler_fix.elf
Propeller Version 1 on com7
Writing 588 bytes to Propeller RAM.
Verifying ... Upload OK!

/* no comment */ ;-)

David Betz · 2011-12-15 07:24

By the way, if you'd like to save some typing you can use -p7 instead of -pcom7 if you want.

Reinhard · 2011-12-16 00:56

David Betz wrote: »

By the way, if you'd like to save some typing you can use -p7 instead of -pcom7 if you want.

Thank you for this tip, because I'm very lazy, I use it.

Reinhard

Dave Hein · 2011-12-16 11:19

Reinhard,

Your program looks very nice. It makes me wonder how much can be packed into a cog C program. It would be interesting to have a contest to see who could produce the largest or most interesting C program that fits in a cog.

Dave

Daniel Harris · 2011-12-22 12:25

We've been considering different contest ideas to encourage testing - especially once we move into beta. I'll add this idea to the list.

Back to the thread, Reinhard, thank you for posting this! It is neat to see complex math being done on the Propeller. Complex being more than just algebra

.

"Graphical" solve of nonlinear differentialequation with the propeller

Comments