% tfyhjennetään ensiksi muisti kaikesta vanhasta roskasta clear all; close all; % ALUSTUS: % tarkastellaan aaltoja suljetussa L-kokoisessa laatikossa, määritellään % koordinaatisto siten, että laatikon vasen reuna on x=-L/2 ja oikea reuna % x=L/2. Jaetaan x-akseli pieniin 0.01-kokoisiin väleihin. % Määrittelemme siis vektorin (listan) xvals, jonka komponentit ovat vain kunkin % paikkaavaruuden pisteen x-koordinaatti, eli % xvals = (-L/2,-L/2+0.1,...,L/2). L=100.0; xvals=-L/2.0:0.01:L/2.0; % yllä muuttuja x on nyt itseasiassa vektori, jonka elementit saavat arvot % -L/2, -L/2+0.01, -L/2+0.02, ..., L/2, % eli siinä on lueteltu kaikki tarkastelemamme x-akselin pisteet % määritellään seuraavaksi gaussinen aaltopakettimme y_x = 1.0/sqrt(2.0*pi)*exp(-0.5*xvals.^2); % myös y_x on nyt vektori, jonka n:s elementti kertoo funktion arvon % vastaavassa x-akselin n:ssä pisteessä. % piirretään vielä kuvaaja figure % luodaan uusi kuvaaja-ikkuna % FOURIER-MUUNNOS % seuraavaksi tehdään aaltopaketille Fourier-muunnos, tai tarkemmin sanoen % cos-muunnos. Koska aaltopaketti on parillinen (sama muoto ja merkki % x:n nollakohdan molemmilla puolilla), katoavat kaikki sin-funktiot % Fourier muunnoksesta. % määritellään ensiksi Fourier muunnoksessa laskemamme n-indeksin arvot. % summaushan menee periaatteessa äärettömyyteen asti, mutta yhä korkeampia % n-indeksejä vastaavat kertoimet \beta_n pienenevät nopeasti. Käytetään % sataa Fourier-termiä (plus vakiotermi, eli n=0, sillä cos(2*pi/L*0*x) = 1). % Määrittelemme siis vektorin nvals, jossa on vain lueteltuna nuo % n-indeksit, eli nvals = (0,1,2,...,100). nvals=0:1:100; % sitten tehdään varsinainen Fourier muunnos, eli jokaista n-indeksin arvoa % nval kohden lasketaan integraali % \beta_nval = (2/L)*\int_{-L/2}^{L/2} f(x)*cos(2*pi/L*nval*x) dx % % integrointi toteutetaan trapz(x,y(x))-funktiolla, joka laskee y(x):n % integraalin yli muuttujan x trapezoidi-menetelmällä. Kyseessä on määrätty % integraali, jossa integroimisalueena on x-muuttujan määrittelyväli. % tehdään silmukka joka juoksee läpi kaikkien n-indeksien arvojen yli beta = zeros(1,length(nvals)); %alustetaan sopivan kokoinen vektori, sisällöllä ei väliä for nindex=1:length(nvals) nval = nvals(nindex); beta(nindex) = (2.0/L)*trapz(xvals,y_x.*cos(2.0*pi/L*nval*xvals)); end % edellä Fourier kertoimet on esitetty indeksin n funktiona. Mutta indeksillä on fysikaalinen merkityksensäkin, % sillä kutakin arvoa n vastaa harmonisen aallon aaltovektori k_n = 2*pi/L*n. Esitetään nyt aaltopaketti % paikka-avaruudessa sekä sen Fourier muunnos aaltovektorin k_n funktiona. kvals = 2.0*pi/L*nvals; % FOURIER-MUUNNOKSEN AIKAKEHITYS % tiedämme miten kukin alkeisaalto cos(pi/L*nval*x) käyttäytyy % ajanfunktiona -- ne kaikki oskilloivat vain ajassa taajuudella, joka % riippuu alkeisaallon aallonpituudesta tai siis n-indeksistä nval % syvän veden aaltojen taajuus riippuu aallonpituudesta näin % f_syva(\lambda) = \sqrt(g/(2\pi))*1/sqrt(\lambda). % matalan veden aaltojen taajuus taasen riippuu näin % f_matala(\lambda) = \sqrt(g*h)/\lambda,% missä g on putoamiskiihtyvyys 9.81 m/s^2 ja h on veden korkeus. Oletetaan % esim. h = 1m % n-indeksiä vastaava aallonpituus puolestaan on \lambda_n = \pi*L/n, % saamme n-indeksiä vastaavan alkeisaallon taajuudeksi gacc = 9.81; %putoamiskiihtyvyys h0 = 1.0; %veden syvyys matalan veden tapauksessa freq_syva_n = sqrt(gacc/(2.0*pi*L)).*sqrt(nvals); freq_matala_n = sqrt(gacc*h0)/(L).*nvals; % ratkaistaan nyt alkuun matalan veden aaltojen leviäminen freq_n = freq_syva_n; %KÄÄNTEISMUUNNOS %lopuksi haluamme tietää miltä aaltopaketti näyttää alkuperäisessä %paikka-avaruudessa, eli x:n funktiona. Meidän pitää tehdä siis käänteinen %Fourier muunnos, tai siis yksinkertaisesti hyödyntää alkuperäistä % funktion esitystä % f(x,t) = \sum_n \beta_n(t) cos(\pi*n/L*x) % y_xt = xvals; %alustetaan sopivan kokoinen vektori, sisällöllä ei väliä title('Aaltopaketin aikakehitys') for t=0.0:0.05:20.0 % aikakehitys beta_t = beta.*cos(2.0*pi*freq_n*t); % käänteismuunnos for xindex=1:length(xvals) xval = xvals(xindex); y_xt(xindex) = sum(beta_t.*cos(2.0*pi/L*nvals*xval)); end % aaltopaketin piirto tiledlayout(2,1) nexttile plot(xvals,y_xt) xlabel('paikka x') ylabel('aallon korkeus y(x)') title(['Aaltopaketti ajanhetkellä t=' num2str(t)]) axis([-L/2 L/2 -0.5 0.5]) nexttile plot(nvals,beta_t,'o') hold on plot(nvals,beta,'.') plot(nvals,-beta,'.') xlabel('n-indeksi') ylabel('Fourier kerroin beta_n(t)') title(['Aaltopaketin Fourier muunnos ajanhetkellä t=' num2str(t)]) hold off pause(0.0625) end