Uživatelské nástroje

Nástroje pro tento web


vyuka:spliny_ve_flexu

Spliny a další interpolace v ActionScriptu

zdrojový soubor ke stažení

package mathExtended
{
	public class MathE
	{
 
 
		public static const NATURAL:Number=9e999;
 
 
 
		//********** pomocná procedura**********************************************
		private static function Locate(xa:Array,x:Number,ilo:int,ihi:int):Object
		{
			 // vstup: 
			 //     xa: pole hodnot
			 //     ilo:  nejnižší index pole xa s hodnotama
			 //     ihi: nejvyšší index pole s hodnotama
			 //     x:  hledaná hodnoty
			 // výstup (pomocí položek objektu) 
			 //    object.khi - dolní index
			 //    object.klo - horní index
			 //    object.h - rozdíl hodnot poležek s horním a dolním indexem
			var outp:Object=new Object(); //návratový objekt
			var k:int;
			//var ihi:int=xa.length-1;
			outp.klo=ilo;
			outp.khi=ihi;
			while (outp.khi-outp.klo>1){
				 k = Math.floor((outp.khi+outp.klo) / 2); //celočíselné dělení
				 if (xa[k] > x){
				 	outp.khi = k; 
				 } else{ 
				 	outp.klo = k; 
				 }		
			}
			outp.h = xa[outp.khi]-xa[outp.klo];
			if (outp.h==0.0){
				trace("Chyba v proceduře SPLINT - špatnně zadané hodnoty vstupních hodnot"); 
			}
			return outp;
		}
		//**********************************************************************************
 
		public static function pchip(x:Array,y:Array,xVal:Number,y1Min:Number=NATURAL,y1Max:Number=NATURAL,iMin:int=0,iMax:int=0):Number
		// funkce je obdobná funkci pchip v Matlabu
		//	vstupy:
		// 		x, y - pole souřadnic uzlových bodů,
		//			pole x a y musí být setříděna podle vzrůstajících hodnot v x 
		//          (pole jsou indexována implicitně od 0, ale parametrem iMin je to možno nastavit od 1)
		// 		xVal - hodnota argumentu, pro který chceme znát interpolovanou hodnotu 
		// 		y1min y1max - požadované derivace v krajních bodech
		// 			pokud y1Min a y1Max se rovnají konstantě NATURAL, jde o přirozený spline
		// 			kdy se hodnoty druhých derivací v krajních bodech se rovnají nule
		//			(implicitně jsou nastaveny hodnoty na NATURAL)
		// 		iMin - index prvního elementu pole, kde začínají hodnoty (implicitně =0)
		//		iMax - index posledního prvku pole, který obsahuje hodnoty - implicitně index posledního prvku pole 
		//				je-li zadáno iMax=0, znamená to totéž jako když iMax = index posledního prvku pole
				// výstup: 
		//    	funkce vrací interpolovanou hodnotu	
		//*****************************************************************	
		//********** vlastní algoritmus výpočtu ***************************
		{
			var y2:Array=Spline(x,y,y1Min,y1Max,iMin,iMax);
		    return SplineInt(x,y,y2,xVal,iMin,iMax);	
		}
 
		public static function Spline(x:Array,y:Array,y1Min:Number=NATURAL,y1Max:Number=NATURAL,iMin:int=0,iMax:int=0):Array
		//	vstupy:
		// 		x, y - pole souřadnic uzlových bodů,
		//			pole x a y musí být setříděna podle vzrůstajících hodnot v x
		// 		y1min y1max - požadované derivace v krajních bodech
		// 			pokud y1Min a y1Max se rovnají konstantě NATURAL, jde o přirozený spline
		// 			kdy se hodnoty druhých derivací v krajních bodech se rovnají nule
		//			(implicitně jsou nastaveny hodnoty na NATURAL)
		// 		iMin - index prvního elementu pole, kde jsou hodnoty (implicitně =0)	
		//		iMax - index posledního prvku pole, který obsahuje hodnoty - implicitně index posledního prvku pole 
		//				je-li zadáno iMax=0, znamená to totéž jako když iMax = index posledního prvku pole
		// výstup: 
		//		pole obsahujícíé druhé derivace v interpolovanmých bodech		
		//*****************************************************************	
		//********** vlastní algoritmus výpočtu ***************************
		{
			if (iMax==0){
				iMax=x.length-1;
			}
		    var k:int;
		    var i:int;
		    var help:Array=new Array();
		    var sig:Number;
		    var p:Number;
		    var qn:Number;
		    var un:Number;
		    var y2:Array=new Array(); //výstupní pole
		    if( iMin==iMax){
                   trace("CHYBA: Nelze provádět interpolaci s jedním bodem");
                   y2=null;
                   return y2;
      		}
 
            if(y1Min >= NATURAL){
 
      			y2[iMin] = 0.0; 
      			help[iMin] = 0.0        	
            }
            else{
      			y2[iMin] = -0.5;
      			help[iMin] = (3.0/(x[iMin+1]-x[iMin]))*((y[iMin+1]-y[iMin])/(x[iMin+1]-x[iMin])-y1Min)       	
            }
            for ( i=iMin+1; i<iMax; i++){
            	sig = (x[i]-x[i-1])/(x[i+1]-x[i-1]);
      			p = sig*y2[i-1]+2.0;
      			y2[i] = (sig-1.0)/p;
      			help[i] = (y[i+1]-y[i])/(x[i+1]-x[i])-(y[i]-y[i-1])/(x[i]-x[i-1]);
      			help[i] = (6.0*help[i]/(x[i+1]-x[i-1])-sig*help[i-1])/p
            }    
			if (y1Max >= NATURAL) {
	  			qn = 0.0; 
	  			un = 0.0			
			}
			else{
	 			qn = 0.5;
      			un = (3.0/(x[iMax]-x[iMax-1]))*(y1Max-(y[iMax]-y[iMax-1])/(x[iMax]-x[iMax-1]));			
			}
   			y2[iMax] = (un-qn*help[iMax-1])/(qn*y2[iMax-1]+1);
   			for( k = iMax-1;k>= iMin;k--){
      			y2[k] = y2[k]*y2[k+1]+help[k];
      		}
 
      		return y2;
		}
		//*************************************************************************************
 
	  	 public static function SplineDerivace(xa:Array,y:Array,y2:Array,xVal:Number,iMin:int=0,iMax:int=0):Number
		// vatupy:
	  	// 		xa, y - pole souřadnic uzlových bodů,
		// 			pole x a y musí být setříděna podle vzrůstajících hodnot v x
		// 		y2 pole obsahujícíé druhé derivace v interpolovanmých bodech - výsledek funkce Spline
		// 		xVal - hodnota argumentu, pro který chceme znát derivaci interpolované hodnoty 
		// 		iMin - index prvního elementu pole, kde jsou hodnoty (implicitně =0)	
		//		iMax - index posledního prvku pole, který obsahuje hodnoty - implicitně index posledního prvku pole 
		//				je-li zadáno iMax=0, znamená to totéž jako když iMax = index posledního prvku pole
		// výstup:
		//    	funkce vrací derivaci v bodě interpolované hodnoty
		//*****************************************************************************************
		//********** vlastní algoritmus výpočtu ***************************************************
	  	 {
			var xx:Number;
			var a:Number;
			var b:Number;
			var o:Object;
			if (iMax==0){
				iMax=xa.length-1;
			}
		    o=Locate(xa,xVal,iMin,iMax);// vrací object.khi,object.klo,object.h)
		   	a = (xa[o.khi]-xVal)/o.h;
	   		b = (xVal-xa[o.klo])/o.h;
	   		if(b<0) { 
	   			b=0;
	   			a=1;
	   		}
	   		else if (a<0) {
	   			a=0;
	   			b=1;
	   		}
	        return (y[o.khi]-y[o.klo])/o.h +((3*b*b-1)*y2[o.khi]*o.h-(3*a*a-1)*y2[o.klo]*o.h)/6.0;
		}
 
 
 
	  	 public static function SplineInt(x:Array,y:Array,y2:Array,xVal:Number,iMin:int=0,iMax:int=0):Number
	  	 {
		// vatupy:
	  	// 		x, y - pole souřadnic uzlových bodů,
		// 			pole x a y musí být setříděna podle vzrůstajících hodnot v x
		// 		y2 pole obsahujícíé druhé derivace v interpolovanmých bodech - výsledek funkce Spline
		// 		xVal - hodnota argumentu, pro který chceme znát interpolovanou hodnotu 
		// 		iMin - index prvního elementu pole, kde jsou hodnoty (implicitně =0)	
		//		iMax - index posledního prvku pole, který obsahuje hodnoty - implicitně index posledního prvku pole 
		//				je-li zadáno iMax=0, znamená to totéž jako když iMax = index posledního prvku pole
		// výstup:
		//    	funkce vrací interpolovanou hodnotu
		//*****************************************************************************************
		//********** vlastní algoritmus výpočtu ***************************************************
			var klo:int;
			var khi:int;
			var xx:Number
			var a:Number;
			var b:Number;
			if (iMax==0){
				iMax=x.length-1;
			}
			var o:Object;
	        o=Locate(x,xVal,iMin,iMax);// vrací object.khi,object.klo,object.h)
	        a = (x[o.khi]-xVal)/o.h;
	        b = (xVal-x[o.klo])/o.h;
	        if (b<0) {
	           return y[iMin]+SplineDerivace(x,y,y2,x[iMin],iMin,iMax)*(xVal-x[iMin]);
	        }
	        else if (a<0) {
	           return y[iMax]+SplineDerivace(x,y,y2,x[iMax],iMin,iMax)*(xVal-x[iMax]);
	        }
	        else {
	           return a*y[o.klo]+b*y[o.khi]+((a*a*a-a)*y2[o.klo]+(b*b*b-b)*y2[o.khi])*(o.h*o.h)/6.0;
	        }
		}
 
 
 
		 public static function SplineRexx(A:Array,B:Array,DY1:Number,DYN:Number):Array
	     {
	     	//upraveno dle Vladimíra Zábranského: Album algoritmů a technik
	     	//http://www.geocities.com/siliconvalley/garage/3323/cz_aat/a_spli.html#spline
	     	// pracuje obdobně jako Spline, dává stejný výsledek
	     	// POZOR POLE ČÍSLOVÁNO OD 1,tj. první hodnoty jsou v A[1] a B[1]
	    	var B2:Array=new Array();
	    	var U:Array=new Array();
	    	var N:int=A.length-1;
	    	var J:int;
	    	var Jm1:int;
	    	var Jp1:int;
	    	var Sig:Number;
	    	var P:Number;
	    	var Nm1:int;
	    	var Qn:Number;
	    	var Un:Number;
	    	if (DY1>0.99E+99){
	    		B2[1]=0;
	    		U[1]=0;
	    	}
	    	else
	    	{
	    		B2[1]=-0.5;
	    	    U[1]=(3/(A[2]-A[1]))*((B[2] - B[1]) / (A[2] - A[1]) - DY1); 
	    	}
	    	for (J=2; J<=N-1;J++){
	    	    Jm1 = J - 1; 
	    	    Jp1 = J + 1;
  				Sig = (A[J] - A[Jm1]) / (A[Jp1] - A[Jm1]);
  				P = Sig * B2[Jm1] + 2; 
  				B2[J] = (Sig - 1) / P; 
  				U[J] = (B[Jp1] - B[J]) / (A[Jp1] - A[J]) -  (B[J] - B[Jm1]) / (A[J] - A[Jm1]);                                                    
  				U[J] = (6 * U[J] / (A[Jp1] - A[Jm1]) - Sig * U[Jm1]) / P 
	    	}
 
	    	Nm1 = N - 1; 
			if (DYN > 0.99E+99) {  
				Qn = 0; 
				Un = 0; 
			}
  			else 
  			{ 
    			Qn = 0.5; 
    			Un = (3 / (A[N] - A[Nm1])) * (DYN - (B[N] - B[Nm1]) / (A[N] - A[Nm1])) 
  			}
 
			B2[N] = (Un - Qn * U[Nm1]) / (Qn * B2[Nm1] + 1);
			for (J = N - 1; J>=1 ;J--) {
  				Jp1 = J + 1 ;
  				B2[J] = B2[J] * B2[Jp1] + U[J] ;
 
			}
			return B2;
	     } 
 
 
	     public static function SplintRexx(A:Array,B:Array,B2:Array,V:Number):Number
	     {
	     	//Toto je také dle Vladimíra Zábranského: Album algoritmů a technik
	     	//http://www.geocities.com/siliconvalley/garage/3323/cz_aat/a_spli.html#spline
	     	// pracuje obdobně jako SplineInt, ale chybně (chybně je naprogramován o v Rexxu!)
	     	//POZOR CHYBNÝ!!!!!
	     	// Až bude čas tak to já nebo někdo další opraví i v Rexxu - ten Zábranského album je velmi 
	     	// užitečný a stojí zato ho opravit (asi původně tam byl správný algoritmus, protože Zábranský uvádí
	     	// správné hodnoty testovacího příkladu - ale algoritmus procedury Splint je implementován chybně) 
	     	var N:int=A.length-1;
	     	var L:int;
	     	var H:int;
	     	var K:int;
	     	var X:Number;
	     	var BV:Number;
	     	var Q:Number;
	     	var Y:Number;
 
	     	L = 1; 
	     	H = N; 
			while ((H - L) > 1){ 
			  K = Math.floor((H + L) / 2); 
			  if (A[K] > V) {
			  	 H = K;
			  }
			  else{
			    L = K;
			  }  
			}
 
			Q = A[H] - A[L];
 
			if (Q == 0) { 
              trace("SPLINT: Vstupní posloupnost musí být rostoucí!!!");
              return NaN;	
  			 }
 
			X = (A[H] - V) / Q; 
			Y = (V - A[L]) / Q; 
 
			BV = X * B[L] + Y * B[H];
 
			BV = BV + (Q*Q) * (X*X*X - X) * B2[L] + (Y*Y*Y - Y) * B2[H] / 6; 
			return BV;   			  			 		     
		}     
	}
 
 
 
}	
vyuka/spliny_ve_flexu.txt · Poslední úprava: 2008/07/24 00:31 autor: kofranek