Jump to content
Sign in to follow this  
aeson

More Math Problems

Recommended Posts

I've tried so many different combinations but I can't seem to get this simple iterated equation to work reliably. It all works fine in visual C++ but things go crazy on the PIC (16F877A).

 

I've actually run into a couple problems, and though I've" worked around" the first problem, perhaps if I get true solution to it, it may help with the second problems.

 

I'm working with a partial differential equations involving the change in temperature along a thickness in material over time (everything is in signed long int):

 

Tnew = T + F0*(T[i-1] - 2*T + T[i+1])/MTPL2;

 

Where MTPL is just a number I multiplied through by and in this case divided through since floats aren't available yet.

 

Putting this in gives me an immediate bad answer. Someone suggested that I also try this to account for arithmatic noise due to integer division truncation:

 

Tnew = T + (F0*T[i-1])/MTPL2 - (2*F0*T)/MTPL2 + (F0*T[i+1])/MTPL2

 

but that didn't work either

 

Just for shits and giggle I even just tried the first two terms:

 

Tnew = T + (F0*T[i-1])/MTPL

 

which still resulted in a huge number as if somewhere the interger looped around. There are a couple other combinations that I used but they didn't work either.

 

That's problem number one. I ended up working around this by putting each term into it's own temporary long int

 

i.e.:

 

temp1 = (F0*T[i-1])/MTPL2;

temp2 = (2*F0*T)/MTPL2;

temp3 = (F0*T[i+1])/MTPL2;

Tnew = T + temp1 - temp2 +temp3;

 

And that seemed to work...for a number of iterations. Why does that work and not the other straight forward method? I thought the order of math operations was established inherently already so this should happen naturally.

 

Still after a number of iterations I get a funky output from out of the blue sometimes it's a huge number, sometimes its a negative number...all I know is that it's not supposed to be there. The funky value appears after the same number of iterations unless I change the code a little, but the code I change usually has nothing to do with the equation, but rather the output to rs232 in an attempt to debug the code. Can anyone help? Here's a snippet of the code below. I can send the whole thing if that will help solve the problem.

 

 

		for(k=6; k<151;k++)
	{
		for(i=1; i<10; i++)
		{
			temp = (F0*T[i-1])/MTPL2;
			temp1 = (2*F0*T[i])/MTPL2;
			temp2 = (F0*T[i+1])/MTPL2;
			Tnew[i] = T[i] + temp - temp1 + temp2;

			Send_Dec(i);
			Send_String(")  ");	
			Send_Dec(T[i]);
			Send_String("  ");		
			Send_Dec(temp);
			Send_String("  ");
			Send_Dec(temp1);
			Send_String("  ");
			Send_Dec(temp2);
			Send_String("  ");
			Send_Dec(T[i] + temp);
			Send_String("  ");
			Send_Dec(T[i] + temp - temp1);
			Send_String("  ");
			Send_Dec(T[i] + temp - temp1 + temp2);
			Send_String("  ");	
			Send_Dec(Tnew[i]);
			Send_String("\n\r");
		}
		Send_String("\n\r");

		Tnew[10] = Tnew[8];

		Send_Dec(k);
		Send_String(")  ");	
		Send_Dec(T[0]/1000);
		Send_String(" ");
		Send_Dec(T[1]/1000);
		Send_String(" ");
		Send_Dec(T[2]/1000);
		Send_String(" ");
		Send_Dec(T[3]/1000);
		Send_String(" ");
		Send_Dec(T[4]/1000);
		Send_String(" ");
		Send_Dec(T[5]/1000);
		Send_String(" ");
		Send_Dec(T[6]/1000);
		Send_String(" ");
		Send_Dec(T[7]/1000);
		Send_String(" ");
		Send_Dec(T[8]/1000);
		Send_String(" ");
		Send_Dec(T[9]/1000);
		Send_String(" ");
		Send_Dec(T[10]/1000);
		Send_String("\n\r\n\r");



		if(Tnew[9] >= 42*MTPL)//bank1 float Tmax = 42; //threshold temp
			{		
				flag = 1;
			}

		for(i= 1; i< 11; i++)
			{
				T[i] = Tnew[i];
			}
	}

Share this post


Link to post
Share on other sites

aeson,

 

can you provide a complete segment of code, ie one that can be compiled and that generates a 'bad' result.

 

Regards

Dave

Share this post


Link to post
Share on other sites
aeson,

 

can you provide a complete segment of code, ie one that can be compiled and that generates a 'bad' result.

 

Regards

Dave

Sure, here it is. I'll just smack the .h file stuff at the top.

When I ran this one, the numbers went funky on iteration 56 the first time I ran it, then at iteration 24, the second and third time I ran it. Thank you so much for any help.

 

Sorry, about all the comments, they help me keep track of everything.

 

 

#include <system.h>	
#include <adc.h>



void serial_Init(char rateScaler);
void Send_Char(char value);
void Send_String(const char* text);
void Send_Dec(signed long n);

signed long T[11] = {0};

signed long Told[11]={0};
signed long Tnew[11] = {0};
//float T0 = 37.0;  //called in program
//bank1 float Ta = 21;//outside ambient temperature during rest
//bank1 float Tmax = 42; //threshold temp
signed long x;
//float A = 5E-08;//this should be between 4*E-8 and 6*E-8
//float L = 0.01;//this is the thickness of the gear in meters
//float S = 15.0; //this is fixed
//float Dt = 0.5; //the time step is fixed
//float F0 = A * Dt * (S/L)*(S/L);
signed long MTPL = 1000000;
signed long MTPL2 = 1000;
signed long F0 =  50; //5E-08 * 0.5 * 1500.0*1500.0;*MTPL
signed long temp, temp1, temp2, temp3;

int i;
int k;
int inc;
int esp;
int flag = 0;

 #pragma DATA 	_CONFIG,  _WDT_OFF & _HS_OSC & _LVP_OFF & _CP_OFF & _PWRTE_OFF





void pause(void){

for(i=0; i<4000; i++);
//for(i=0; i<4000; i++);
//for(i=0; i<4000; i++);
//for(i=0; i<4000; i++);


};

void init_a2d(void)
{
trisa = 11111111b;
porta = 0;
trisc = 0;
portc = 0;
trisb = 0;
portb = 0;
//	CMCON0 = 7; //turn off comparators in 16F684
adcon0=0;	// select Fosc/2
adcon1=0;	// select left justify result. A/D port configuration 0
set_bit(adcon0 , ADON);		// turn on the A2D conversion module //	ANSEL = 1; in 16F684
clear_bit( intcon, GIE );
};

/* Return an 8 bit result */

unsigned char read_a2d(void){
	int tmp;
set_bit(adcon0 , ADON);	
//ADCON0=0b00000001;	// clear current channel select
set_bit(adcon0, GO_DONE);	
while(adc_go)continue;
return(adresh);	// return 8 MSB of the result
};


void main(void)
{

init_a2d();
   serial_Init(64);//9600
portb = 0b01010101;
Send_String("\fStartUp...\n\r");
T[0] = 27*MTPL;
T[1] = 28*MTPL;
T[2] = 29*MTPL;
T[3] = 30*MTPL;
T[4] = 31*MTPL;
T[5] = 32*MTPL;
T[6] = 33*MTPL;
T[7] = 34*MTPL;
T[8] = 35*MTPL;
T[9] = 36*MTPL;
T[10] = 35*MTPL;

  	while(1)
  	{
	portc = 0b11111111;

	x = read_a2d();//adc_measure(1);

	T[0] = MTPL*x/2;
	Send_String("Ts=");             
	Send_Dec(T[0]);
	Send_String("\n\r");
	for(k=1;k<6;k++)
	{

		for(i=1;i<10;i++)
		{

			temp = (F0*T[i-1])/MTPL2;
			temp1 = (2*F0*T[i])/MTPL2;
			temp2 = (F0*T[i+1])/MTPL2;
			Tnew[i] = T[i] + temp - temp1 + temp2;
			Send_String("  ");
/*
			temp = -2*T[i];
			temp = temp+T[i+1];
			temp = T[i-1]+temp;
			temp = temp/MTPL2;
			temp = F0*temp;
			Tnew[i] = T[i]+temp;
*/
		}
		Send_String("\n\r");

		Tnew[10] = Tnew[8];

/*
		i = 1;
		Send_Dec(T[i]);
		Send_String("  ");		
		Send_Dec(temp);
		Send_String("  ");
		Send_Dec(temp1);
		Send_String("  ");
		Send_Dec(temp2);
		Send_String("  ");
		Send_Dec(T[i] + temp - temp1 + temp2);
		Send_String("\n\r\n\r");

		Send_Dec(T[i]);
		Send_String("  ");		
		Send_Dec((F0*T[0])/MTPL2);
		Send_String("  ");
		Send_Dec((2*F0*T[1])/MTPL2);
		Send_String("  ");
		Send_Dec((F0*T[2])/MTPL2);
		Send_String("  ");
		Send_Dec(T[1] + (F0*T[0])/MTPL2 - (2*F0*T[1])/MTPL2 + (F0*T[2])/MTPL2);
		Send_String("\n\r\n\r");
*/

		for(i = 1; i<11; i++)
			T[i] = Tnew[i];

		Send_Dec(T[0]/1000);
		Send_String("  ");
		Send_Dec(T[1]/1000);
		Send_String("  ");
		Send_Dec(T[2]/1000);
		Send_String("  ");
		Send_Dec(T[3]/1000);
		Send_String("  ");
		Send_Dec(T[4]/1000);
		Send_String("  ");
		Send_Dec(T[5]/1000);
		Send_String("  ");
		Send_Dec(T[6]/1000);
		Send_String("  ");
		Send_Dec(T[7]/1000);
		Send_String("  ");
		Send_Dec(T[8]/1000);
		Send_String("  ");
		Send_Dec(T[9]/1000);
		Send_String("  ");
		Send_Dec(T[10]/1000);
		Send_String("\n\r\n\r");



	}

	//Store Temperatures
	for(k = 1; k<11;k++)
		Told[k] = T[k];

	T[0] = 21*MTPL;//outside ambient temperature during rest

	//Elapse 150*Dt/60 minutes for prediction
	for(k=6; k<151;k++)
	{
		for(i=1; i<10; i++)
		{
			temp = (F0*T[i-1])/MTPL2;
			temp1 = (2*F0*T[i])/MTPL2;
			temp2 = (F0*T[i+1])/MTPL2;
			Tnew[i] = T[i] + temp - temp1 + temp2;

			Send_Dec(i);
			Send_String(")  ");	
			Send_Dec(T[i]);
			Send_String("  ");		
			Send_Dec(temp);
			Send_String("  ");
			Send_Dec(temp1);
			Send_String("  ");
			Send_Dec(temp2);
			Send_String("  ");
			Send_Dec(T[i] + temp);
			Send_String("  ");
			Send_Dec(T[i] + temp - temp1);
			Send_String("  ");
			Send_Dec(T[i] + temp - temp1 + temp2);
			Send_String("  ");	
			Send_Dec(Tnew[i]);
			Send_String("\n\r");
		}
		Send_String("\n\r");

		Tnew[10] = Tnew[8];

		Send_Dec(k);
		Send_String(")  ");	
		Send_Dec(T[0]/1000);
		Send_String(" ");
		Send_Dec(T[1]/1000);
		Send_String(" ");
		Send_Dec(T[2]/1000);
		Send_String(" ");
		Send_Dec(T[3]/1000);
		Send_String(" ");
		Send_Dec(T[4]/1000);
		Send_String(" ");
		Send_Dec(T[5]/1000);
		Send_String(" ");
		Send_Dec(T[6]/1000);
		Send_String(" ");
		Send_Dec(T[7]/1000);
		Send_String(" ");
		Send_Dec(T[8]/1000);
		Send_String(" ");
		Send_Dec(T[9]/1000);
		Send_String(" ");
		Send_Dec(T[10]/1000);
		Send_String("\n\r\n\r");



		if(Tnew[9] >= 42*MTPL)//bank1 float Tmax = 42; //threshold temp
			{		
				flag = 1;
			}

		for(i= 1; i< 11; i++)
			{
				T[i] = Tnew[i];
			}
	}

	for(k = 1; k < 11; k++)
		T[k] = Told[k];

}
}

void serial_Init(char rateScaler)
{
   spbrg = rateScaler;                                                 //baud rate init
   txsta = 0b00100100;                                                  //transmit init
   rcsta = 0b10010000;                                                  //receive init
}

void Send_Char(char value)
{
   while((txsta & 1 << TRMT) == 0);                                 // TRMT is better then TXIF
   txreg = value;
}

void Send_String(const char* text)
{
   char i = 0;
   while(text[i] != 0)
       Send_Char(text[i++]);
}

void Send_Dec(signed long n)
{
if(n<0)
{
	n = n*(-1);
	Send_String("-");
}
   if(n > 999999999)
       Send_Char(((n / 1000000000) % 10) + '0');
   if(n > 99999999)
       Send_Char(((n / 100000000) % 10) + '0');
   if(n > 9999999)
       Send_Char(((n / 10000000) % 10) + '0');
   if(n > 999999)
       Send_Char(((n / 1000000) % 10) + '0');
   if(n > 99999)
       Send_Char(((n / 100000) % 10) + '0');
   if(n > 9999)
       Send_Char(((n / 10000) % 10) + '0');
   if(n > 999)
       Send_Char(((n / 1000) % 10) + '0');
   if(n > 99)
       Send_Char(((n / 100) % 10) + '0');
   if(n > 9)
       Send_Char(((n / 10) % 10) + '0');
   Send_Char((n % 10) + '0');



}

Share this post


Link to post
Share on other sites

aeson,

 

You provide a little bit too much.

 

Without having your setup its a bit hard to test.

 

Can you just provide the expression with some numbers assigned to the terms that demonstrate the problem, ie so that the code can be run and the problem demonstrated without all the hardware.

 

Something we can run under SourceBoost simulator is what we really want.

 

Regards

Dave

Share this post


Link to post
Share on other sites
aeson,

 

You provide a little bit too much.

 

Without having your setup its a bit hard to test.

 

Can you just provide the expression with some numbers assigned to the terms that demonstrate the problem, ie so that the code can be run and the problem demonstrated without all the hardware.

 

Something we can run under SourceBoost simulator is what we really want.

 

Regards

Dave

 

sorry about that, here it is again with all the fat cut out. I didn't know if or how to watch different indecies of an array so I put everything in "temp1" through "temp10" from the corresponding indexed T's. When I tried to run the simulator myself, the watch list only let my look at T[0] dubbed just plain "T". Everything should progress down to 21000000 after a very long time, and it does for a short time, then everything goes wacky. As I mentioned I have two problems, the other problem is demonstrated by temp11, 12, and 13, note the comments.

 

Thanks again for the help.

 

-Jason

#include <system.h>	
#pragma DATA 	_CONFIG,  _WDT_OFF & _HS_OSC & _LVP_OFF & _CP_OFF & _PWRTE_OFF


signed long T[11] = {0};
signed long Told[11]={0};
signed long Tnew[11] = {0};
signed long MTPL2 = 1000;
signed long F0 =  50; 
signed long tmp, tmp1, tmp2;
signed long temp1, temp2, temp3, temp4, temp5, temp6, temp7;
signed long temp8, temp9, temp10, temp11, temp12, temp13, temp14;

int i;
int k;


void main(void)
{

T[0] = 27000000;
T[1] = 28000000;
T[2] = 29000000;
T[3] = 30000000;
T[4] = 31000000;
T[5] = 32000000;
T[6] = 33000000;
T[7] = 34000000;
T[8] = 35000000;
T[9] = 36000000;
T[10] = 35000000;


//Weird Answers Part 1:
//this is just an example of math that doesn't seem to work
//Going for 
// Answer = 28000000 + 50*27000000/1000
//        = 28000000 + 50*27000
//		  = 28000000 + 1350000
//		  = 29350000

//this works
temp11 = T[1] + F0*T[0]/MTPL2;

//but this doesn't
i = 1;
temp12 = T[i] + F0*T[i-1]/MTPL2;

//but this does...
tmp = F0*T[i-1]/MTPL2;
temp13 = T[i] + tmp;

//why?


//Weird Answers Part 2:
//Works for awhile but then odd seemingly random answers pop up
  	while(1)
  	{


	T[0] = 21000000;//outside ambient temperature during rest


	for(k=1; k<150;k++)
	{
		for(i=1; i<10; i++)
		{

			tmp = (F0*T[i-1])/MTPL2;
			tmp1 = (2*F0*T[i])/MTPL2;
			tmp2 = (F0*T[i+1])/MTPL2;
			Tnew[i]	= T[i] + tmp - tmp1 + tmp2;

			//this is what I'm going for
			//Tnew[i] = T[i]+F0*(T[i-1]-2*T[i]+T[i+1])/MTPL2;


/*
			//something else I had tried
			tmp = -2*T[i];
			tmp = tmp+T[i+1];
			tmp = T[i-1]+tmp;
			tmp = tmp/MTPL2;
			tmp = F0*tmp;
			Tnew[i] = T[i]+tmp;
*/


		}
		Tnew[10] = Tnew[8];


		temp1 = Tnew[1];
		temp2 = Tnew[2];
		temp3 = Tnew[3];
		temp4 = Tnew[4];
		temp5 = Tnew[5];
		temp6 = Tnew[6];
		temp7 = Tnew[7];
		temp8 = Tnew[8];
		temp9 = Tnew[9];
		temp10 = Tnew[10];

		for(i= 1; i< 11; i++)
			{
				T[i] = Tnew[i];
			}
	}

	for(k = 1; k < 11; k++)
		T[k] = Told[k];

}
}

 

 

 

Also, While I'm here posting, I was also wondering why I get an error saying

 

Failed to locate file 'E:\Prg\picant\compiler\picantc\libc\div.c'

Failed to locate file 'E:\Prg\picant\compiler\picantc\libc\mul.c'

 

when I don't use the snapshot option. I don't even have a harddrive e: and I couldn't find a place to change the directory.

Share this post


Link to post
Share on other sites

Looking at your original equation:aeson Posted on: Apr 28 2006, 07:50 PM

Tnew = T + F0*(T[i-1] - 2*T + T[i+1])

 

Where F0 is some constant < 1 and T are points along the line.

This function would be easier to evaluate with integer math by using a few transforms:

 

Eliminate multiply by 2:

 

Tnew = T + F0*(T[i-1] - T + T[i+1] - T)

 

Set F0 as a ratio of two integers: Fn/Fd.

 

Where Fn < Fd, and where Fd*ABS(2*Tmax) <= INTEGER_MAX and where Fn*ABS(4*Tmax) <= INTEGER_MAX.

 

Tnew = (Fd*T + Fn*(T[i-1] - T + T[i+1] - T)) / Fd

 

Doing this should help you identify where and when an arithmetic overflows occurs.

 

Controlling the range of the input T values is essential to the calculation stability.

Values outside of a fairly narrow range will cause an arithmetic overflow.

 

The limits suggested for the Fn and Fd values assure that arithmetic overflow should not occur for T values within a reasonable range.

 

A good choice for Fn and Fd is: Fn=1, Fd=256

 

This give a Tmax value of 4,194,303.

 

Assuming your input is in degrees C then bias the input values by 10,000 and input range runs from +41.94303C to -41.94303C. Input values outside of this range will cause an arithmetic overflow.

 

If you need execute this calculation quickly you will need to code this in assembler as the C code and 32-bit math function will be inefficient when the Fn and Fd values are not optimized for byte shifts and multiplies and divides are used.

Share this post


Link to post
Share on other sites

aeson,

 

Looks like a compiler bugs :)

 

This:

Tnew = T + tmp - tmp1 + tmp2;

does not evalauted correctly.

 

The longer hand version:

Tnew = T;

Tnew += tmp;

Tnew -= tmp1;

Tnew +=tmp2;

does evaute correctly.

 

Also looking at your other problem, which is similarly related, its seems that more that one array index in an expression with 32bit maths causes the problem.

 

Regards

Dave

Share this post


Link to post
Share on other sites

Here is a version of your code that seems to work correctly with integer math:

#include <system.h>

// Target 16F877A

// Configuration word

#pragma DATA     _CONFIG,  _WDT_OFF & _HS_OSC & _LVP_OFF & _CP_OFF & _PWRTE_OFF

#pragma CLOCK_FREQ 20000000

/*

Demo program to iterate the partial differential equation:

  T[i]new = T[i] + F0*(T[i-1] - 2*T[i] + T[i+1])
  
Temporary variables are used to work around a compiler bug
that hits calculation using 32-bit integer arrays.

Using temporary variables does have the added benefit that
the PDE can be iterated in place.

TODO: 

Run time check interim calculations for arithmetic overflow.

*/

signed long T[11];

#define BIAS (10000)
#define Fn (1)
#define Fd (16)
#define ITERATE_LIMIT (150)

void main(void)
{
signed long tmp0, tmp1, tmp2, absErr;
unsigned char i;
unsigned short k;

// load test data
T[0] = 27*BIAS;
T[1] = 28*BIAS;
T[2] = 29*BIAS;
T[3] = 30*BIAS;
T[4] = 31*BIAS;
T[5] = 32*BIAS;
T[6] = 33*BIAS;
T[7] = 34*BIAS;
T[8] = 35*BIAS;
T[9] = 36*BIAS;
T[10] = 35*BIAS;

while(1)
{
  T[0] = 21*BIAS;  //outside ambient temperature during rest
  k=0;
  do
  {
    //T[i] = (Fd*T[i]+Fn*(T[i-1]-T[i]+T[i+1]-T[i]))/Fd
    absErr = 0;
    tmp1 = T[0];
    tmp2 = T[1];
    for(i=1; i<10; i++)
    {
      // set tmp0 = T[i-1], tmp1 = T[i], tmp2 = T[i+1]
      tmp0 = tmp1;
      tmp1 = tmp2;
      tmp2 = T[i+1];

      // compute error term Fn*(T[i-1]-T[i]+T[i+1]-T[i])
      tmp0 -= tmp1;
      tmp0 += tmp2;
      tmp0 -= tmp1;
      tmp0 *= Fn;

      // update array T[i] = (Fd*T[i]+Fn*(T[i-1]-T[i]+T[i+1]-T[i]))/Fd

      // this does not work
      //T[i] *= Fd;
      //T[i] += tmp0;
      //T[i] /= Fd;

      // this works
      T[i] = (T[i] * Fd + tmp0) / Fd;

      // sum absolute values of all error terms
      if (tmp0 < 0) absErr -= tmp0;
      else absErr += tmp0;
    }
    
    T[10] = T[8]; // I am not sure why you are doing this

    k++; // count iteration

    // exit early if we have iterated to a small enough error.
    if (absErr < 2) break;

  } while ((k != 0) && (k < ITERATE_LIMIT));

  nop(); // just for a debug breakpoint

 }
}

Edited by cac001

Share this post


Link to post
Share on other sites

Join the conversation

You are posting as a guest. If you have an account, sign in now to post with your account.
Note: Your post will require moderator approval before it will be visible.

Guest
Reply to this topic...

×   Pasted as rich text.   Paste as plain text instead

  Only 75 emoji are allowed.

×   Your link has been automatically embedded.   Display as a link instead

×   Your previous content has been restored.   Clear editor

×   You cannot paste images directly. Upload or insert images from URL.

Loading...
Sign in to follow this  

×
×
  • Create New...