Uogólniony (symetryczny) problem własny - wyznaczanie modów własnych struny w 1D

Uogólniony (symetryczny) problem własny - wyznaczanie modów własnych struny w 1D

Tomasz Chwiej 26 marca 2020

1 Wprowadzenie

Na laboratorium zajmiemy się wyznaczaniem częstości drgań własnych struny, której wychylenie w czasie i przestrzeni opisuje funkcja ψ = ψ(x, t). Dynamiką struny rządzi równanie falowe (N-naciąg struny, ρ(x)- liniowy rozkład gęstości):

N ρ(x)

∂²ψ

∂x² = ∂²ψ

∂t² (1)

Dokonujemy separacji zmiennych: i) najpierw podstawiając ψ(x, t) = u(x)θ(t), a następnie ii) dzieląc przez iloczyn uθ

N ρ(x)

1 u

∂²u

∂x² = 1 θ

∂²θ

∂t² = const =−λ (λ = ω², ω− częstość własna drgań) (2) dzięki czemu otrzymujemy równanie różniczkowe zależne tylko od zmiennej położeniowej

−∂²u

∂x² = λρ(x)

N u (3)

Struna przymocowana jest w punktach ±L/2 (L-długość struny). Wprowadzamy siatkę równoodle- głych węzłów: x = xi, u(x) = ui, ρ(x) = ρi Odległość pomiędzy węzłami wynosi

∆x = L

n + 1 (4)

a położenie w przestrzeni wyznaczamy tak xi =−L

2 + ∆x· (i + 1), i = 0, 1, 2, . . . , n − 1 (5) Teraz możemy dokonać dyskretyzacji równania (3) podstawiając trójpunktowy iloraz różnicowy cen- tralny za drugą pochodną

−ui−1− 2ui+ ui+1

∆x² = λρi

Nui (6)

co można zapisać w postaci (A,B - macierze, uuu - wektor)

Auuu = λBuuu (7)

co stanowi tzw.uogólniony problem własny, w którym elementy macierzowe są zdefiniowane następu- jąco

Ai,j = (−δi,j+1+ 2δi,j− δi,j−1) /∆x² (8)

1

oraz

Bi,j = ρ_i

Nδi,j (9)

gdzie

δ_i,j =

{ 1, i = j

0, i̸= j (10)

jest deltą Kroneckera.

2 Zadania do wykonania

1. Przyjmujemy następujące parametry: L = 10, n = 200, ρ(x) = 1 + 4 α x², N = 1.

2. Utworzyć macierze A i B oraz wypełnić je zgodnie z wzorami (8) i (9).

3. Rozwiązać równanie (7) dla α ∈ [0, 100] z krokiem ∆α = 2. Dla każdej wartości parametru α do pliku zapisać wartości pierwiastków z 6 kolejnych najmniejszych wartości własnych i sporządzić odpowiedni wykres (ω =√

λ = f (α)).

4. Dla α = 0 oraz α = 100 zapisać do pliku wektory własne odpowiadające 6 najniższym wartościom własnym i sporządzić ich wykresy.

Przykładowe wyniki

0 0.05 0.1 0.15 0.2 0.25

0 20 40 60 80 100

ω

α (a) ω₁ ω₂ ω₃ ω₄ ω₅ ω₆

-0.1 -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08 0.1

-4 -2 0 2 4

u(x)

x (b) α=0

u₁ u₂ u₃ u₄ u₅ u₆

-0.15 -0.1 -0.05 0 0.05 0.1 0.15

-4 -2 0 2 4

u(x)

x (c) α=100

u₁ u₂ u₃ u₄ u₅ u₆

Rysunek 1: (a) zmiany wartości własnych w funkcji parametru α (uwaga: zakres osi pionowej ograni- czono do 0.25) oraz wykresy wektorów własnych odpowiadających 6 najniższym wartościom własnym dla α = 0 (b) i α = 100 (c).

Dla α = 0 (rys. 1(b)) wektory własne są nieparzystymi wielokrotnościami ”połówek” sinusa. Na wykresie 1(c) widać natomiast że dla α = 100, środek struny staje się masywny przez co wektory 1,3,5 mają płaską centralną część.

3 Uwagi

1. W projekcie korzystamy z biblioteki GSL dołączając pliki nagłówkowe

#include</usr/include/gsl/gsl_eigen.h>

2. Dla każdej wartości α macierze zawsze wypełniamy zgodnie z wzorami (8) i (9), łącznie z zerami (eliminujemy w ten sposób śmieci numeryczne z poporzednich przebiegów pętli).

3. Macierze A i B są symetryczne, więc do rozwiązania uogólnionego problemu własnego używamy metody GSL-a

int gsl_eigen_gensymmv(gsl_matrix * A, gsl_matrix * B, gsl_vector * eval, gsl_matrix * evec, gsl_eigen_gensymm_workspace * w) gdzie: eval to wektor wartości własnych (nieposortowany), evec to macierz n× n w której ko- lumnach zapisane są wektory własne, a wektor pomocniczy w definiujemy tak

gsl_eigen_gensymmv_workspace *w = gsl_eigen_gensymmv_alloc(n);

4. Wartości i wektory własne sortujemy (po rozwiązaniu problemu) stosując funkcję GSL-a int gsl_eigen_gensymmv_sort(gsl_vector * eval, gsl_matrix * evec,

gsl_eigen_sort_t sort_type);

gdzie podstawiamy

gsl eigen sort t=GSL EIGEN SORT ABS ASC (sortowanie od najmniejszej do największej wartości własnej)