#08 ; Lagrange-interpolation

Nu skall vi börja interpolera på riktigt. Vi kommer att använda de tabeller för tiden t, sträckan s som kulan flugit och kulans hastighet v som vi har behandlat förr.

Lagrange-interpolation från 4 värden är ganska enkel att implementera för argument vars avstånd från varandra är inte det samma. Vi kommer att interpolera från sträckor i tabeller som x och hastigheter eller tider som y. Tider i tabellen är ju i jämna mellanrum, men sträckor och hastigheter är det inte. Principen är enkel, vi kommer att presentera relationen y = f(x) som finns gömd i tabeller, som en tredje grads polynom med några konstanter a₀, a₁, a₂ samt a₃ vars värde är tills vidare okänt.

y = a0· x3 + a1· x2 + a2· x + a3

Den sträcka för vilken vi vill ha de andra värden heter x och resultaten heter y. Det är emellertid bättre att presentera formeln först i en "cubisk" form. De fyra argument (sträcka) från tabellen kommer att heta x₀, x₁, x₂ och x₃ och de andra värden (hastigheten eller tiden) kommer att heta y₀, y₁, y₂ och y₃.

y = A · ( x - x1 ) · ( x - x2 ) · ( x - x3 )
  + B · ( x - x0 ) · ( x - x2 ) · ( x - x3 )
  + C · ( x - x0 ) · ( x - x1 ) · ( x - x3 )
  + D · ( x - x0 ) · ( x - x1 ) · ( x - x2 )

Men hur mycket är A, B, C och D? Man kan lösa tex. A genom att ersätta x = x₀ och y = y₀ i termen där det finns A

y0 = A · ( x0 - x1 ) · ( x0 - x2 ) · ( x0 - x3 )
A = y0 / ( ( x0 - x1 ) · ( x0 - x2 ) · ( x0 - x3 ) )

Och likadant för B, C och D. Den färdiga formeln set ut så här:

y = y0 · ( x - x1 )·( x - x2 )·( x - x3 ) / ( ( x0 - x1 )·( x0 - x2 )·( x0 - x3 ) )
  + y1 · ( x - x0 )·( x - x2 )·( x - x3 ) / ( ( x1 - x0 )·( x1 - x2 )·( x1 - x3 ) )
  + y2 · ( x - x0 )·( x - x1 )·( x - x3 ) / ( ( x2 - x0 )·( x2 - x1 )·( x2 - x3 ) )
  + y3 · ( x - x0 )·( x - x1 )·( x - x2 ) / ( ( x3 - x0 )·( x3 - x1 )·( x3 - x2 ) )

Formeln skulle bli mycket enklare om all x-värden våre på samma avstånd från varandra, men det är nu inte fallet. Vi får nöja oss med den där långa formeln. Det kan se ut litet förvirrande i början, men man hittar där en viss logik när man granskar det närmare. Vi behöver ju inte minnas formeln utantill eller använda det på egen hand. Programmet skall göra det arbetsamma jobbet för oss. Oss räcker att vi förstår vad som pågår.

Samma princip kunde man utnyttja med längre tabeller, men det kan vi i praktiken inte rekommendera. En Lagrange-polynom av hög grad kan mycket lätt producera missvisande resultat, på grund av små fel i basinformationen. Metoden tvingar polynomen gå exakt genom givna punkter och då kan hända att en metod av hög grad hittar en resultat som visserligen är rätt rent matematiskt, men i praktiken en mycket konstig och orealistisk lösning. En polynom av för hög grad kan "vibrera" om det finnst litet felaktiga siffror i tabeller. I praktiken kan de värden från vilka man skall interpolera, inte vara helt exakta.

Nästa steget i vår kära lilla projekt kommer att bli det att jag implementerar Lagrange-algoritmen i JavaScript och förenar det med den numerisk integration som producerar kulans hastighet och sträckan. Först skall jag tänka hur det borde låta som pseudokod. Det blir något annorlunda än jag tänkte förr, därför att vi interpolerar från 4 och inte från 3 värden.

// Funktionen kallas två gånger per en viss sträcka (önskadeX), 
// först med tabeller s[] och v[] ; vi får hastigheten
// och sedan med tabeller s[] och t[] ; vi får tiden

FUNKTION Lagrange3djeGrad ( xtabel[], ytabel[], önskadeX )

  // Ack vad dumt, men det hjälper inte

  x_x0 := önskadeX - xtabel[0]
  x_x1 := önskadeX - xtabel[1]
  x_x2 := önskadeX - xtabel[2]
  x_x3 := önskadeX - xtabel[3]

  x0_x1 := xtabel[0] - xtabel[1]
  x0_x2 := xtabel[0] - xtabel[2]
  x0_x3 := xtabel[0] - xtabel[3]
  x1_x0 := -x0_x1 
  x2_x0 := -x0_x2
  x3_x0 := -x0_x3

  x1_x2 := xtabel[1] - xtabel[2]
  x1_x3 := xtabel[1] - xtabel[3]
  x2_x1 := -x1_x2
  x3_x1 := -x1_x3

  x2_x3 := xtabel[2] - xtabel[3]
  x3_x2 := -x2_x3

  y0 := ytabel[0] · x_x1 · x_x2 · x_x3 / ( x0_x1 · x0_x2 · x0_x3 )
  y1 := ytabel[1] · x_x0 · x_x2 · x_x3 / ( x1_x0 · x1_x2 · x1_x3 )
  y2 := ytabel[2] · x_x0 · x_x1 · x_x3 / ( x2_x0 · x2_x1 · x2_x3 )
  y3 := ytabel[3] · x_x0 · x_x1 · x_x2 / ( x3_x0 · x3_x1 · x3_x2 )

  interpoleradY := y0 + y1 + y2 + y3

  RETURN ( interpoleradY )
SLUT FUNKTION

Då mäste vi prova på det, här är programmet.

Meny
Huvudsida (finska)