Znajdowanie wielomianowych aproksymacji fali sinusoidalnej

Chcę aproksymować sinusoidę podaną przez poprzez zastosowanie wielomianowego wavehapera do prostej fali trójkątnej , generowanej przez funkcję$\sin\left(\pi x\right)$

$$T\left(x\right)=1-4\left|\tfrac{1}{2}-\operatorname{mod}(\tfrac{1}{2}x+\tfrac{1}{4},\ 1)\right|$$

gdzie nazwa jest ułamkową częścią :$\operatorname{mod}(x, 1)$$x$

$$\operatorname{mod}(x, y) \triangleq y \cdot \left( \left\lfloor \frac{x}{y}\right\rfloor - \frac{x}{y} \right)$$

Seria Taylor może być wykorzystana jako zmieniacz fal.

$$S_1\left(x\right)=\frac{\pi x}{2}-\frac{\frac{\pi x}{2}^3}{3!}+\frac{\frac{\pi x}{2}^5}{5!}-\frac{\frac{\pi x}{2}^7}{7!}$$

Biorąc pod uwagę powyższe funkcje, $S_1(T(x))$ da nam przyzwoite przybliżenie fali sinusoidalnej. Ale musimy przejść do siódmej potęgi serii, aby uzyskać rozsądny wynik, a szczyty są nieco niskie i nie będą miały nachylenia dokładnie zero.

Zamiast serii Taylora, moglibyśmy użyć wielomianu zmieniającego fale zgodnie z kilkoma zasadami.

• Musi przejść przez -1, -1 i + 1, + 1.
• Nachylenie przy -1, -1 i + 1, + 1 musi wynosić zero.
• Musi być symetryczny.

Prosta funkcja spełniająca nasze wymagania:

$$S_2\left(x\right)=\frac{3x}{2}-\frac{x^3}{2}$$

Wykresy $S_2(T(x))$ i $\sin\left(\pi x\right)$ są dość bliskie, ale nie tak bliskie jak seria Taylora. Pomiędzy szczytami i przejściami przez zero wyraźnie się nieco odbiegają. Cięższa i dokładniejsza funkcja spełniająca nasze wymagania:

$$S_3\left(x\right)=\frac{x(x^2-5)^2}{16}$$

Prawdopodobnie jest to wystarczająco blisko dla moich celów, ale zastanawiam się, czy istnieje inna funkcja, która bardziej zbliża się do fali sinusoidalnej i jest tańsza obliczeniowo. Bardzo dobrze rozumiem, jak znaleźć funkcje spełniające trzy powyższe wymagania, ale nie jestem pewien, jak znaleźć funkcje spełniające te wymagania, a także najlepiej pasujące do fali sinusoidalnej.

Jakie istnieją metody znajdowania wielomianów, które naśladują falę sinusoidalną (po zastosowaniu do fali trójkątnej)?

---

Aby wyjaśnić, niekoniecznie szukam tylko nieparzystych wielomianów, chociaż są to najprostszy wybór.

Coś w rodzaju następującej funkcji może również pasować do moich potrzeb:

$$S_4\left(x\right)=\frac{3x}{2}+\frac{x^2}{4}+\frac{x^4}{4}$$

Spełnia to wymagania w zakresie ujemnym i można zastosować rozwiązanie fragmentaryczne, aby zastosować go również w zakresie dodatnim; na przykład

$$\frac{3x}{2}-\frac{P\left(x,2\right)}{4}-\frac{P\left(x,4\right)}{4}$$

gdzie jest podpisaną funkcją mocy .$P$

Byłbym również zainteresowany rozwiązaniami wykorzystującymi podpisaną funkcję mocy do obsługi wykładników ułamkowych, ponieważ daje to nam kolejne „pokrętło do skręcania” bez dodawania kolejnego współczynnika.

$$a_0x\ +a_1P\left(x,\ p_1\right)$$

Biorąc pod uwagę odpowiednie stałe, mogłoby to potencjalnie uzyskać bardzo dobrą dokładność bez ciężkości wielomianów piątego lub siódmego rzędu. Oto przykład spełniający opisane tutaj wymagania przy użyciu wybranych ręcznie: .$a_0=1.\overline{666}, a_1=-0.\overline{666}, p_1=2.5$

$$\frac{5x-2P\left(x,\ \frac{5}{2}\right)}{3}$$

W rzeczywistości te stałe są bardzo zbliżone do i , i . Podłączenie ich daje coś, co wygląda wyjątkowo blisko fali sinusoidalnej.$\frac \pi 2$$1 - \frac \pi 2$$e$

$$\frac{\pi}{2}x\ +\left(1-\frac{\pi}{2}\right)P\left(x,e\right)$$

Innymi słowy, wygląda bardzo blisko między 0,0 a π / 2,1. Wszelkie przemyślenia na temat tego znaczenia? Może narzędzie takie jak Octave może pomóc odkryć „najlepsze” stałe dla tego podejścia.$x-\frac{x^e}{6}$$\sin(x)$

około dziesięć lat temu zrobiłem to dla nienazwanej firmy zajmującej się syntezatorami muzyki, która miała badania i rozwój niedaleko mojego mieszkania w Waltham MA. (nie mogę sobie wyobrazić, kim oni są.) Nie mam współczynników.

ale spróbuj tego:

$$\begin{align} f(x) &\approx \sin\left(\tfrac{\pi}{2} x \right) \qquad \qquad \text{for }-1 \le x \le +1 \\ \\ &= \tfrac{\pi}{2} x (a_0 + a_1 x^2 + a_2 x^4) \\ \end{align}$$

gwarantuje to, że .$f(-x) = -f(x)$

Aby zagwarantować, że wtedy$f'(x) \Big|_{x=\pm1} =0$

$$f'(x) = \tfrac{\pi}{2}( a_0 + 3 a_1 x^2 + 5 a_2 x^4)$$

$$a_0 + 3 a_1 + 5 a_2 = 0 \tag{1}$$

To pierwsze ograniczenie. Aby zagwarantować, że zatem$|f(\pm 1)| = 1$

$$a_0 + a_1 + a_2 = \tfrac{2}{\pi} \tag{2}$$

To drugie ograniczenie. Eliminowanie i rozwiązywanie Eqs. (1) i (2) dla w kategoriach (które należy dostosować):a 2 a $a_0$$a_2$$a_1$

$$a_0 = \tfrac{5}{2 \pi} - \tfrac{1}{2} a_1$$

$$a_2 = -\tfrac{1}{2 \pi} - \tfrac{1}{2} a_1$$

Teraz masz tylko jeden współczynnik, , który możesz , aby uzyskać najlepszą wydajność:$a_1$

$$f(x) = \tfrac{\pi}{2} x \left((\tfrac{5}{2 \pi} - \tfrac{1}{2} a_1) + a_1 x^2 - (\tfrac{1}{2 \pi} + \tfrac{1}{2} a_1)x^4 \right)$$

W ten sposób przekręciłbym aby uzyskać najlepszą wydajność oscylatora fali sinusoidalnej. Zmieniłbym użycie powyższego i symetrii fali sinusoidalnej około i umieściłem dokładnie jeden cały cykl w buforze o mocy dwóch liczb punktów (powiedzmy 128, nie obchodzi mnie to) i uruchomiłem FFT na tym idealny cykl. x = $a_1$$x=1$

Bin 1 wyniku FFT będzie siłą sinusoidy i powinien wynosić około . Teraz możesz wyregulować aby zwiększyć lub zmniejszyć zniekształcenie trzeciej harmonicznej. od , aby . To jest w bin 3 wyników FFT. Jednak zniekształcenie piątej harmonicznej (wartość w bin 5) będzie następowało (będzie rosło wraz ze spadkiem trzeciej harmonicznej). Dostosowałbym aby siła 5. poziomu harmonicznej była równa 3. poziomowi harmonicznej. Będzie to około -70 dB od pierwszej harmonicznej (jak pamiętam). To będzie najlepiej brzmiąca fala sinusoidalna z taniego, 3-współczynnikowego, nieparzysto-symetrycznego wielomianu 5. rzędu.a 1 a 1 ≈ $N/2$$a_1$a0≈1a$a_1 \approx \tfrac{5}{\pi} - 2$$a_0 \approx 1$$a_1$

Ktoś inny może napisać kod MATLAB. Jak ci to brzmi?

Zwykle dokonuje się aproksymacji minimalizującej pewną normę błędu, często L_ -norm (gdzie minimalny błąd maksymalny jest minimalizowany) lub norma (gdzie minimalny średni błąd kwadratu). - aproksymacja odbywa się za pomocą algorytmu wymiany Remeza . Jestem pewien, że możesz znaleźć kod Open Source implementujący ten algorytm. Jednak w tym przypadku uważam, że bardzo prosta (dyskretna) l_2 jest wystarczająca. Spójrzmy na kod Matlab / Octave i wyniki: L 2 L ∞ l $L_{\infty}$$L_2$$L_{\infty}$$l_2$

x = przestrzeń linowa (0, pi / 2300); % siatki w [0, pi / 2]
x = x (:);
% przekroczonego układu równań liniowych
% (tylko przy użyciu mocy nieparzystych)
A3 = [x, x. ^ 3];
A5 = [x, x. ^ 3, x. ^ 5];
b = sin (x);
% rozwiązania w sensie l2
c3 = A3 \ b;
c5 = A5 \ b;
f3 = A3 * c3; % Przybliżenie trzeciego rzędu
f5 = A5 * c5; % Przybliżenie 5. rzędu

Poniższy rysunek pokazuje błędy aproksymacji dla zamówienia i dla aproksymacji . Błędy maksymalne zbliżanie się i odpowiednio. 5 t $3^{rd}$$5^{th}$8.8869e-031.5519e-04

Optymalne współczynniki to

c3 =
   0,988720369237930
  -0,149993929056091

i

c5 =
   0,99976918199047515
  -0,16582163562776930
   0,00757183954143367

Przybliżenie trzeciego rzędu to

$$\sin(x)\approx 0.988720369237930x-0.144993929056091x^3,\quad x\in[-\pi/2,\pi/2]\tag{1}$$

a przybliżeniem piątego rzędu jest

$$\sin(x)\approx 0.99976918199047515x-0.16582163562776930x^3+\\0.00757183954143367x^5,\quad x\in[-\pi/2,\pi/2]\tag{2}$$

EDYTOWAĆ:

Spojrzałem na przybliżenia z podpisaną funkcją mocy, jak zasugerowano w pytaniu, ale najlepsze przybliżenie nie jest prawie lepsze niż przybliżenie trzeciego rzędu pokazane powyżej. Funkcja aproksymująca to

$$f(x)=x - \frac{1}{p}\left(\frac{\pi}{2}\right)^{1-p}x^p,\qquad x\in[0,\pi/2]\tag{3}$$

gdzie stałe zostały wybrane w taki sposób, i . Moc została zoptymalizowana, aby osiągnąć najmniejszy błąd maksymalny w zakresie . Stwierdzono, że optymalna wartość wynosi . Poniższy rysunek pokazuje błędy aproksymacji dla aproksymacji trzeciego rzędu i dla nowej aproksymacji :f ' ( π / 2 ) = 0 p [ 0 , π / 2 ] p p = 2,774 ( 1 ) ( 3 $f'(0)=1$$f'(\pi/2)=0$$p$$[0,\pi/2]$$p$$p=2.774$$(1)$$(3)$

Maksymalny błąd aproksymacji aproksymacji jest , ale należy pamiętać, że aproksymacja trzeciego rzędu przekracza ten błąd tylko w pobliżu i że w większości jego błąd aproksymacji jest w rzeczywistości mniejszy niż błąd podpisanej funkcji mocy .π /$(3)$4.5e-3$\pi/2$

EDYCJA 2:

Jeśli nie przeszkadza ci podział, możesz również użyć wzoru przybliżenia sinusoidy Bhaskara I , który ma maksymalny błąd przybliżenia 1.6e-3:

$$\sin(x)\approx\frac{16x(\pi-x)}{5\pi^2-4x(\pi-x)},\qquad x\in [0,\pi/2]\tag{4}$$

Rozpocznij od sparametryzowanego wielomianu nieparzystej symetrii 5. rzędu:

$$\begin{align} f(x) &= a_0 x^1 + a_1 x^3 + a_2 x^5 \\ &= x\big( a_0 + a_1 x^2 + a_2 x^4 \big) \\ &= x\bigg( a_0 + x^2 \Big(a_1 + a_2 x^2\Big) \bigg) \\ \end{align}$$

Teraz nakładamy pewne ograniczenia na tę funkcję. Amplituda powinna wynosić 1 na szczytach, innymi słowy $f(1) = 1$ . Podstawienie $1$ na $x$ daje:

$$a_0 + a_1 + a_2 = 1 \tag1$$

To jedno ograniczenie. Nachylenie na szczytach powinno wynosić zero, innymi słowy $f'(1) = 0$ . Pochodna $f(x)$ to

$$a_0 + 3 a_1 x^2 + 5 a_2 x^4$$

i podstawienie $1$ na $x$ daje nasze drugie ograniczenie:

$$a_0 + 3 a_1 + 5 a_2 = 0 \tag2$$

Teraz możemy wykorzystać nasze dwa ograniczenia do rozwiązania za $a_1$ i $a_2$ w kategoriach $a_0$ .

$$a_1=\frac{5}{2}-2a_0 \\ a_2=a_0-\frac{3}{2}\tag{3}$$

Pozostaje to, aby dostosować $a_0$ aby uzyskać ładny dopasowanie. Nawiasem mówiąc, 0 (i nachylenie przy początku) kończy się ≈ gatunku$a_0$$\approx\frac{\pi}{2}$ , jak widać zwykresufunkcji.

Optymalizacja parametrów

Poniżej przedstawiono szereg optymalizacji współczynników, które skutkują tymi względnymi amplitudami harmonicznych w porównaniu do częstotliwości podstawowej (1. harmonicznej):

W złożonej serii Fouriera :

$$\sum_{k=-\infty}^\infty c_k e^{i\tfrac{2\pi}{P} kx},$$

rzeczywistego $P\text{-}$ okresowego przebiegu o $P = 4$ i symetrii czasowej około $x = 1$ i przy połowie okresu zdefiniowanego przez funkcję nieparzystą $f(x)$ powyżej $-1 \le x \le 1,$ współczynnik $k\text{th}$ złożonej harmonicznej wykładniczej wynosi:

$$c_k = \frac{1}{P}\int_{-1}^{-1+P}\left(\cases{f(x)&if $x < 1$\\-f(x-2)&if $x \ge 1$}\right)e^{-i\tfrac{2\pi}{P}kx}dx.$$

Z powodu zależności $2 \cos(x) = e^{ix} + e^{-ix}$ (patrz: wzór Eulera ), amplituda rzeczywistej harmonicznej sinusoidalnej przy $k > 0$ wynosi $2\left|c_k\right|,$ która jest dwa razy większa niż złożoność wykładnicza tej samej częstotliwości. Można to wmasować do postaci, która ułatwia symbolicznemu oprogramowaniu matematycznemu uproszczenie całki:

$$2|c_k| = \frac{2}{4}\left|\int_{-1}^{3}\left(\cases{f(x)&if $x < 1$\\-f(x-2)&if $x \ge 1$}\right)e^{-i\tfrac{2\pi}{4}kx}dx\right|\\ = \frac{1}{2}\left|\int_{-1}^{1}f(x)e^{-i\tfrac{\pi}{2}kx}dx - \int_{1}^{3}f(x-2)e^{-i\tfrac{\pi}{2}kx}dx\right|\\ = \frac{1}{2}\left|\int_{-1}^{1}f(x)e^{-i\tfrac{\pi}{2}kx}dx - \int_{-1}^{1}f(x+2-2)e^{-i\tfrac{\pi}{2}k(x+2)}dx\right|\\ = \frac{1}{2}\left|\int_{-1}^{1}f(x)e^{-i\tfrac{\pi}{2}kx}dx - \int_{-1}^{1}f(x)e^{-i\tfrac{\pi}{2}k(x+2)}dx\right|\\ = \frac{1}{2}\left|\int_{-1}^{1}f(x)\left(e^{-i\tfrac{\pi}{2}kx} - e^{-i\tfrac{\pi}{2}k(x+2)}\right)dx\right|\\ = \frac{1}{2}\left|e^{i\tfrac{\pi}{2}x}\int_{-1}^{1}f(x)\left(e^{-i\tfrac{\pi}{2}kx} - e^{-i\tfrac{\pi}{2}k(x+2)}\right)dx\right|\\ = \frac{1}{2}\left|\int_{-1}^{1}f(x)\left(e^{-i\frac{\pi}{2}k(x-1)}-e^{-i\frac{\pi}{2}k(x+1)}\right)dx\right|$$

Powyższe korzysta z tego $|e^{ix}|=1$ dla prawdziwego $x.$Niektóre systemy algebry komputerowej łatwiej upraszczają całkę, zakładając, że $k$ jest rzeczywiste, i upraszczają liczbę całkowitą $k$ na końcu. Wolfram Alpha może całkować poszczególne warunki całki końcowej odpowiadające warunkom wielomianu $f(x)$ . Dla współczynników podanych w równaniu. 3 otrzymujemy amplitudę:

$$= \left|\frac{48\left((-1)^k - 1\right)\left(16\,a_0\left(\pi^2 k^2 - 10\right) - 5\times(5\pi^2 k^2 - 48)\right)}{\pi^6k^6}\right|$$

5. rząd, pochodna ciągła

Możemy rozwiązać za wartość $a_0$ zapewniający równe amplitudzie $2|c_k|$trzeciej i piątej harmonicznej. Będą dwa rozwiązania odpowiadające trzeciej i piątej harmonicznej o równych lub przeciwnych fazach. Najlepszym rozwiązaniem jest to, które minimalizuje maksymalną amplitudę trzeciej i wyższej harmonicznej oraz równoważnie maksymalną względną amplitudę trzeciej i wyższej harmonicznej w porównaniu do częstotliwości podstawowej (1. harmonicznej):

$$a_0 = \frac{3\times(132375\pi^2 - 130832)}{16\times(15885\pi^2 - 16354)}\approx 1.569778813,\\ a_1 = \frac{5}{2} - 2a_0 = \frac{79425\pi^2 - 65416}{8\times(-15885\pi^2 + 16354)}\approx -0.6395576276,\\ a_2 = a_0 - \frac{3}{2} = \frac{15885\pi^2}{16\times(15885\pi^2 - 16354)}\approx 0.06977881382.$$

Daje to częstotliwość podstawową przy amplitudzie $\frac{13679616}{15885\pi^6 - 16354\pi^4} \approx 1.000071420$oraz zarówno 3., jak i 5. harmonicznej przy amplitudzie względnej$\frac{1}{8906}$ lub około$-78.99\text{ dB}$porównaniu do częstotliwości podstawowej. $k\text{th}$harmoniczna ma względna amplituda$\frac{\left(1 - (-1)^k\right)\left|8177k^2 - 79425\right|}{142496 k^6}.$

7. rząd, pochodna ciągła

Podobnie optymalne przybliżenie wielomianowe siódmego rzędu z tymi samymi ograniczeniami początkowymi oraz harmoniczną 3., 5. i 7. na najniższym możliwym równym poziomie wynosi:

$$\begin{align} f(x) &= a_0 x^1 + a_1 x^3 + a_2 x^5 + a_3 x^7\\ &= x\big( a_0 + a_1 x^2 + a_2 x^4 + a_3 x^7 \big) \\ &= x\bigg( a_0 + x^2 \Big(a_1 + x^2 \big(a_2 + a_3 x^2 \big) \Big) \bigg) \\ \end{align}$$

$$a_0 = \frac{2a_2 + 4a_3 + 3}{2} \approx 1.570781972,\\ a_1 = -\frac{4a_2 + 6a_3 + 1}{2} \approx -0.6458482979,\\ a_2 = \frac{347960025\pi^4 - 405395408\pi^2}{16\times(281681925 \pi^4 - 405395408 \pi^2 + 108019280)} \approx 0.07935067784,\\ a_3 = - \frac{16569525\pi^4}{16\times(281681925\pi^4 - 405395408\pi^2 + 108019280)} \approx -0.004284352588.$$

Jest to najlepsze z czterech możliwych rozwiązań odpowiadających kombinacjom równych / przeciwnych faz trzeciej, piątej i siódmej harmonicznej. Częstotliwość podstawowa ma amplitudę $\frac{2293523251200}{281681925\pi^8 - 405395408\pi^6 + 108019280\pi^4} \approx 0.9999983752,$a trzecia, piąta i siódma harmoniczna mają amplitudę względną$\frac{1}{1555395} \approx -123.8368\text{ dB}$porównaniu do podstawowego. $k\text{th}$harmoniczna ma względna amplituda$\frac{\left(1 - (-1)^k\right)\left|1350241k^4 - 50674426k^2 + 347960025\right|}{597271680k^8}$ porównaniu do podstawowego.

5. zamówienie

Jeśli wymóg ciągłej pochodnej zostanie porzucony, przybliżenie 5. rzędu będzie trudniejsze do symbolicznego rozwiązania, ponieważ amplituda 9. harmonicznej wzrośnie powyżej amplitudy 3., 5. i 7. harmonicznej, jeśli są one ograniczone równe i zminimalizowane. Testując 16 różnych rozwiązań odpowiadających różnym podzbiorom trzech harmonicznych z $\{3, 5, 7, 9\}$ o jednakowej amplitudzie i jednakowych lub przeciwnych fazach, najlepszym rozwiązaniem jest:

$$f(x) = a_0 x^1 + a_1 x^3 + a_2 x^5\\ a_0 = 1 - a_1 - a_2 \approx 1.570034357\\ a_1 = \frac{3\times(2436304\pi^2 - 2172825\pi^4)}{8\times(1303695\pi^4 - 1827228\pi^2 + 537160)} \approx -0.6425216143\\ a_2 = \frac{1303695\pi^4}{16\times(1303695\pi^4 - 1827228\pi^2 + 537160)} \approx 0.07248725712$$

Częstotliwość podstawowa ma amplitudę $\frac{1080430592}{1303695\pi^6 - 1827228\pi^4 + 537160\pi^2} \approx 0.9997773320.$3., 5. i 9. harmoniczna mają amplitudę względną$\frac{7}{263777} \approx -91.52\text{ dB},$a 7. harmoniczna ma amplitudę względną$\frac{726083}{31033100273} \approx -92.6\text{ dB}$porównaniu do podstawowego. $k\text{th}$harmoniczna ma względna amplituda$\frac{\left(1 - (-1)^k\right)\left|67145k^4 - 2740842k^2 + 19555425\right|}{33763456k^6}.$

To przybliżenie ma niewielki róg na granicach półcyklu , ponieważ wielomian ma zerową pochodną nie przy $x = \pm 1$ ale przy $x \approx \pm 1.002039940.$Przy $x = 1$ wartość pochodnej wynosi około $0.004905799828$ . Powoduje to wolniejsze asymptotyczne zanikanie amplitud harmonicznych przy dużym $k,$ porównaniu do aproksymacji 5. rzędu, która ma pochodną ciągłą.

7. zamówienie

Przybliżenie siódmego rzędu bez pochodnej ciągłej można znaleźć podobnie. Podejście to wymaga przetestowania 120 różnych rozwiązań i zostało zautomatyzowane przez skrypt Python na końcu tej odpowiedzi. Najlepszym rozwiązaniem jest:

$$f(x) = a_0 x^1 + a_1 x^3 + a_2 x^5 + a_3 x^7\\ a_0 = 1 - a_1 - a_2 - a_3 \approx 1.5707953785726114835\\ a_1 = - \frac{5\times(4374085272375\pi^6 - 6856418226992\pi^4 + 2139059216768\pi^2)}{16\times(2124555703725\pi^6 - 3428209113496\pi^4 + 1336912010480\pi^2 - 155807094720)} \approx -0.64590724797262922190\\ a_2 = \frac{2624451163425\pi^6 - 3428209113496\pi^4}{16\times(2124555703725\pi^6 - 3428209113496\pi^4 + 1336912010480\pi^2 - 155807094720)} \approx 0.079473610232926783079\\ a_3 = -\frac{124973864925\pi^6}{16\times(2124555703725\pi^6 - 3428209113496\pi^4 + 1336912010480\pi^2 - 155807094720)} \approx -0.0043617408329090447344$$

Częstotliwość podstawowa ma amplitudę $\frac{16991801282396160}{2124555703725\pi^8 - 3428209113496\pi^6 + 1336912010480\pi^4 - 155807094720\pi^2} \approx 1.0000024810802368487.$Największa względna amplituda harmonicznych powyżej podstawy wynosi$\frac{50}{2400688077}\approx -133.627\text{ dB}.$w porównaniu do podstawowego. $k\text{th}$harmoniczna ma względna amplituda$\frac{\left(1 - (-1)^k\right)\left|-162299057k^6 + 16711400131k^4 - 428526139187*k^2 + 2624451163425\right|}{4424948250624k^8}.$

Źródło Python

from sympy import symbols, pi, solve, factor, binomial

numEq = 3 # Number of equations
numHarmonics = 6 # Number of harmonics to evaluate

a1, a2, a3, k = symbols("a1, a2, a3, k")
coefficients = [a1, a2, a3]
harmonicRelativeAmplitude = (2*pi**4*a1*k**4*(pi**2*k**2-12)+4*pi**2*a2*k**2*(pi**4*k**4-60*pi**2*k**2+480)+6*a3*(pi**6*k**6-140*pi**4*k**4+6720*pi**2*k**2-53760)+pi**6*k**6)*(1-(-1)**k)/(2*k**8*(2*pi**4*a1*(pi**2-12)+4*pi**2*a2*(pi**4-60*pi**2+480)+6*a3*(pi**6-140*pi**4+6720*pi**2-53760)+pi**6))

harmonicRelativeAmplitudes = []
for i in range(0, numHarmonics) :
    harmonicRelativeAmplitudes.append(harmonicRelativeAmplitude.subs(k, 3 + 2*i))

numCandidateEqs = 2**numHarmonics
numSignCombinations = 2**numEq
useHarmonics = range(numEq + 1)

bestSolution = []
bestRelativeAmplitude = 1
bestUnevaluatedRelativeAmplitude = 1
numSolutions = binomial(numHarmonics, numEq + 1)*2**numEq
solutionIndex = 0

for i in range(0, numCandidateEqs) :
    temp = i
    candidateNumHarmonics = 0
    j = 0
    while (temp) :
        if (temp & 1) :
            if candidateNumHarmonics < numEq + 1 :
                useHarmonics[candidateNumHarmonics] = j
            candidateNumHarmonics += 1
        temp >>= 1
        j += 1
    if (candidateNumHarmonics == numEq + 1) :
        for j in range(0,  numSignCombinations) :
            eqs = []
            temp = j
            for n in range(0, numEq) :
                if temp & 1 :
                    eqs.append(harmonicRelativeAmplitudes[useHarmonics[0]] - harmonicRelativeAmplitudes[useHarmonics[1+n]])
                else :
                    eqs.append(harmonicRelativeAmplitudes[useHarmonics[0]] + harmonicRelativeAmplitudes[useHarmonics[1+n]])
                temp >>= 1
            solution = solve(eqs, coefficients, manual=True)
            solutionIndex += 1
            print "Candidate solution %d of %d" % (solutionIndex, numSolutions)
            print solution
            solutionRelativeAmplitude = harmonicRelativeAmplitude
            for n in range(0, numEq) :                
                solutionRelativeAmplitude = solutionRelativeAmplitude.subs(coefficients[n], solution[0][n])
            solutionRelativeAmplitude = factor(solutionRelativeAmplitude)
            print solutionRelativeAmplitude
            solutionWorstRelativeAmplitude = 0
            for n in range(0, numHarmonics) :
                solutionEvaluatedRelativeAmplitude = abs(factor(solutionRelativeAmplitude.subs(k, 3 + 2*n)))
                if (solutionEvaluatedRelativeAmplitude > solutionWorstRelativeAmplitude) :
                    solutionWorstRelativeAmplitude = solutionEvaluatedRelativeAmplitude
            print solutionWorstRelativeAmplitude
            if (solutionWorstRelativeAmplitude < bestRelativeAmplitude) :
                bestRelativeAmplitude = solutionWorstRelativeAmplitude
                bestUnevaluatedRelativeAmplitude = solutionRelativeAmplitude                
                bestSolution = solution
                print "That is a new best solution!"
            print

print "Best Solution is:"
print bestSolution
print bestUnevaluatedRelativeAmplitude
print bestRelativeAmplitude

Czy pytasz o to z powodów teoretycznych lub praktycznego zastosowania?

Zwykle, gdy masz kosztowną funkcję obliczeniową w ograniczonym zakresie, najlepszą odpowiedzią jest zestaw tabel odnośników.

Jednym z podejść jest stosowanie najlepiej dopasowanych parabol:

n = podłoga (x * N + .5);

d = x * N - n;

i = n + N / 2;

y = L_0 + L_1 [i] * d + L_2 [i] * d * d;

Znalezienie paraboli w każdym punkcie, który spełnia wartości d wynoszące -1/2, 0 i 1/2, zamiast używania pochodnych w 0, zapewnia ciągłe przybliżenie. Możesz również przesunąć wartość x zamiast indeksu tablicy, aby poradzić sobie z ujemnymi wartościami x.

Ced

=================================================

Kontynuacja:

Ilość wysiłku i wyniki, które zostały włożone w znalezienie dobrych przybliżeń, są imponujące. Byłem ciekawy, jak wyglądałoby moje nudne i nijakie paraboliczne rozwiązanie. Nic dziwnego, że robi się znacznie lepiej. Oto wyniki:

   Metoda Minimalna Maksymalna Średnia RMS
  -------- -------- -------- -------- --------
     Moc -8,48842 1,99861 -4,1936 5,23002
    OP S_3 -2,14675 0,00000 -1,20299 1,40854
     Bhask -1,334370 1,63176 -0,14367 0,97353
     Stosunek -0,24337 0,22770 -0,00085 0,16244
     rbj 5 -0,06724 0,15519 -0,00672 0,04195
    Olli5C -0,16367 0,20212 0,01003 0,1268
     Olli5 -0,26698 0,00000 -0,15177 0,16402
    Olli7C -0,00213 0,00000 -0,00129 0,00143
     Olli7 -0,00005 0,00328 0,00149 0,00181
    Para16 -0,00921 0,00916 -0,00017 0,00467
    Para32 -0,00104 0,00104 -0,00001 0,00053
    Para64 -0,00012 0,00012 -0,00000 0,00006

$\pi/2$

$x-\frac{x^e}{6}$

Linia rbj 5 jest taka sama jak rozwiązanie Matt L5 c5.

16, 32 i 64 to liczba przedziałów, które mają napady paraboliczne. Oczywiście w pierwszej pochodnej występują niewielkie nieciągłości na każdej granicy przedziału. Wartości funkcji są jednak ciągłe. Zwiększenie liczby interwałów tylko zwiększa wymagania dotyczące pamięci (i czas inicjalizacji), nie zwiększa ilości obliczeń potrzebnych do aproksymacji, która jest mniejsza niż w przypadku innych równań. Wybrałem potęgi dwóch, ponieważ implementacja punktu stałego może uratować podział, używając AND w takich przypadkach. Nie chciałem też, aby liczba była proporcjonalna do próbkowania próbnego.

Uruchomiłem program pytlowy Olli Niemitalo i otrzymałem to jako część wydruku: „Rozwiązanie kandydujące 176 ze 120”. Myślałem, że to dziwne, więc wspominam o tym.

Jeśli ktoś chce, żebym włączył inne równania, daj mi znać w komentarzach.

Oto kod częściowych aproksymacji parabolicznych. Cały program testowy jest zbyt długi, aby opublikować.

# ================================================= ============================
def FillParab (argArray, argPieceCount):

# y = ad ^ 2 + bd + c

# ym = a .25 - b. 5 + c
# y = c
# yp = a .25 + b .5 + c

# c = y
# b = yp - ym
# a = (yp + ym - 2 lata) * 2

# ---- Oblicz tablice wyszukiwania

        theStep = pi * .5 / float (argPieceCount - 1)
        theHalf = theStep * .5

        theL0 = zera (argPieceCount)
        theL1 = zera (argPieceCount)
        theL2 = zera (argPieceCount)

        dla k w zakresie (0, argPieceCount):
         x = float (k) * theStep

         ym = sin (x - theHalf)
         y = sin (x)
         yp = sin (x + theHalf)

         theL0 [k] = y
         theL1 [k] = yp - ym
         theL2 [k] = (yp + ym - 2,0 * y) * 2

# ---- Wypełnij

        theN = len (argArray)

        theFactor = pi * .5 / float (theN - 1)

        dla i w zakresie (0, theN):
         x = float (i) * theFactor

         kx = x / theStep
         k = int (kx + .5)
         d = kx - k

         argArray [i] = theL0 [k] + (theL1 [k] + theL2 [k] * d) * d

# ================================================= ============================

=======================================

Dodatek

$S_3$