\documentclass[a4paper]{article}

\usepackage[latin1]{inputenc}
\usepackage{amsmath}
\usepackage[swedish]{babel}
\usepackage[T1]{fontenc}
\usepackage[dvips]{graphics}
\title{Laboration 6 i 2D1250 Tillämpade Numeriska Metoder II}
\author{Adam Jonsson, Arifali Hirji}
\date{}

\begin{document} 
\maketitle 
\section{Poissons ekvation}

Varmefordelningen i stationart tillstånd for vår modell ges av Poissons ekvation: \begin{equation} \begin{split}
\nabla^2 T(x,y)=f(x,y)
\end{split}\end{equation}dar $f(x,y)=-q_0L^2$ om $0.002<x<0.018, 0<y<0.002$ och 0 annars, samt 'Fouriers lag' vilken ger oss villkor på randen:\begin{equation}-k_c\textbf{n}\cdot T(x,y)=h(T-T\star)\end{equation}
Har ar $q_0=5.86\cdot 10^7$. \\Då vi infor $u(x,y)=k_c(T(Lx,Ly)-T\star)$ får vi med variablebytet \textbf{x'}=L\textbf{x} och att $\frac{d}{d\textbf{x} }=L\frac{d}{d\textbf{x'} }$ ur (1) och (2):\begin{equation}\begin{split}u_x_x+u_y_y&=q_0, 0.1<x<0.9, 0<y<0.2\\- u_y(x,0)&=0\\-u_y(x,0.2)&=\tilde{h}u(x,0.2)\\u_x(0,y)&=\tilde{h}u(0,y)\\-u_x(1,y)&=\tilde{h}u(1,y),\textbf{{x}}\in D \end{split}\end{equation}
dar $\tilde{h}=-\frac{Lh)}{k_C_n}=0.0050, h=100$ (allt i SI). \\Vi diskretiserade problemet genom att forst fixera M och på så satt dela upp definitionsområdet i $(N+2)(M+2)$ punkter. Vi kallar losningen i de olika punkterna  $u_k_,_l=u(k\triangle x,l\triangle x), k=0,1,2,..,M+1,l=0,1,2,..,N+1.$ Har ar $\triangle x=\frac{1}{1+N}$. Vi har alltså de $MN$ st inre punkterna:$u_1_1,u_1_2,..u_1_M,u_2_1,...,u_N_M$. \\Det forenklade randvillkoret ger $u_k_l=0$ då $k=0,N+1$ eller $l=0,M+1$. Vi får också $f_k_,_l=f(k\triangle,l\triangle)=q_0$ for inre punkter då $\frac{k}{1+N}k>0.1,\frac{k}{1+N}k<0.9$ och $\frac{l}{1+M}k<0.1$. $f$'s utseende kan ses i figur nedan. Vi noterar att vi kommer att approximera en diskontinuerlig funktion. \\Sedan anvander vi fempunktsapproximationen av $u_x_x +u_y_y$ ... Genom att infora $\textbf{u}=(u_1_1,u_1_2,u_1_3,u_2_1,u_2_2,u_2_3,..,u_9_3)$ och $MN$-vektorn $\textbf{q} _0$ vilken innehåller vardena på $f_k_,_l$, och ar ordnad på samma satt som $u_k_l$ får vi ekvationen for de inre punkterna:
\begin{equation} \textbf{\emph{A}}\textbf{\emph{u}}=\textbf{\emph{q}} _0 \end{equation}
dar $\textbf{\emph{A}}$ ar $MN\timesMN$ och har komponenterna $A_i_j=-4, i=j, A_i_j=1, |i-j|=1$, eller $|i-j|=M$. Resultaten ses i figur 1. Vi korde aven filen chip1, (vilken loser samma problem) i FEMLAB. Den losningen ses i Fig 2.

\begin{figure}[h]
\begin{center}
\scalebox{0.5}{\includegraphics{forenklade_rv.eps}}
\caption{Forenklade randvillkor}
\end{center}
\end{figure}

\begin{figure}[h]
\begin{center}
\scalebox{0.5}{\includegraphics{nr1sv.eps}}
\caption{FEMLAB}
\end{center}
\end{figure}

Då vi byter ut de forenklade randvillkoren mot de verkliga och approximerar derivatorna på randen enl anvisning har vi alltså: $u_{k,0}=u_{k,1},u_{k,M+1}=u_{k,M}$samt $u_{0,l}=\frac{1}{1+\tilde{h}\triangle}u_{1,l},u_{N+1,l}=\frac{1}{1+\tilde{h}\triangle}u_{N,l}$ Detta påverkar endast diagonalelementen $A_j_j$ enl foljande: $\frac{1}{1+\tilde{h}\triangle}$ adderas for $j=1,2,..,M,j=MN-M,..MN$ och dessutom adderas 1 om $j=mM$ eller $j=(m+1)M, m=1,2,..M$. Då vi loser $\textbf{\emph{A}}\textbf{\emph{u}}=\textbf{\emph{q}}$ med denna modifikation får losningen foljande utseende. 

%\subsection{Losning av problemet med FEMLAB}
%Då vi korde de fardiga filerna i FEMLAB fick vi 

\section{} \subsection{Vågekvationen i 1D - Losning med Leap-Frog-Metoden}\\I den har uppgiften studerar vi vågekvationen:
\begin{equation}
\begin{split}
u_t_t-u_x_x&=0,  0<x<1\\u(0,x)&=u_0(x)\\u_t(0,x)&=v_0(x)
\end{split}
\end{equation}
\\Vi anvander Leap-frog-metoden \begin{equation}u_{j}^{n+1}- 2u_{j}^{n} + u_{j}^{n-1}=\lambda^2(u_{j+1}^{n}-2u_{j}^{n} + u_{j-1}^{n},)  \lambda=\frac{\triangle t}{\triangle x} \end{equation} Detta kan ses som en numerisk approximation med centraldifferenser i tiden av ODE systemet  
\begin{equation}\begin{split}
\textbf{\emph{u}}_t_t=\textbf{\emph{A}}\textbf{\emph{u}}
\end{split}
\end{equation}
dar vi infort vektorn $\textbf{u}=(u(t,\triangle x),u(t,2\triangle x),..,u(t,N\triangle x))$. $\textbf{u}$ innehåller alltså losningarna for de de olika punkterna på x-axeln. Vi har också att $A_i_j=-2/(\triangle x)^2$ om $i=j$ och $1/(\triangle x)^2$ om $|i-j|=1$
\\Då vi loser detta system får vi alltså for varje u-komponent: \begin{equation}\begin{split}\frac{(u^{n+1}-2u^n+u^{n-1})_j}{(\triangle t)^2}=(\textbf{\emph{A}}\textbf{\emph{u}})_j\end{split}\end{equation} \\Vi implementerade sedan metoden i MATLAB. I figur XX ses losningen vid x=0.25 med jamna mellanrum i intervallet t=[0 2]
\begin{figure}[h]
\begin{center}
\scalebox{0.5}{\includegraphics{ogon.eps}}
\caption{Ogonblicksbilder t=linspace(0,2,2/dt)}
\end{center}
\end{figure}

Då vi studerade losningen med hjalp av drawnode kunde vi observera att vågorna reflekteras omvanda vid $x=0,x=1$, dar vågen ar 'inspand' (vi har Dirichlet -villkor).

\subsection{Modellering av trycket i en flojt}Vi anvande sedan en variant av vågekvationen for att modellera trycket i en flojt av langd $L$.
\begin{equation}\begin{split}p_t_t-c_0^2p_x_x+D(x)p=G(t,x),  p(0)=p(L)=0\end{split}\end{equation}
dar $c_0=344m/s$, $D(x)$ beror positivt av storleken på fingerhålet vid $x$ och $G(t,x)$ anvands for att modellera inblåsningen i flojten.
\\Efter variabelbyte enl $u(t,x)=p(Lt/c_0,Lx) $ får vi 
\begin{equation}\begin{split}u_t_t-u_x_x+d(x)p=g(t,x),  u(0)=u(1)=0\end{split}\end{equation}
Precis som i 2.1 kan vi se det som att vi loser ODEn 
\begin{equation}\begin{split}\textbf{\emph{u}}_t_t=\textbf{\emph{A}}\textbf{\emph{u}}-\textbf{\emph{D}}\textbf{\emph{u}}+\textbf{\emph{g}}(t)\end{split}\end{equation}
med centraldifferenser i tiden.
\\Vi borjade med att titta på losningen då $d(x)=0$ och $g(t,x)=e^{-10t}sin(\pi x)$. Fig XX visar den stående våg som den driften (inblåsningen i flojten) ger upphov till. \\
\begin{figure}[h]
\begin{center}
\scalebox{0.5}{\includegraphics{staende.eps}}
\caption{}
\end{center}
\end{figure}

Vi studerade sedan hur trycket varierade i punkten $x=0.25$ for $0<t<20$. Den oscillerande losningen visas i Fig XX.

\begin{figure}[h]
\begin{center}
\scalebox{0.5}{\includegraphics{x025.eps}}
\caption{}
\end{center}
\end{figure}

 Vi undersokte sedan sammansattningen hos den stående vågen. Fourieranalysen sager ju att en svangning kan beskrivas som en summa av fundamentala svangningar, dar koefficienterna i serien ar svangningens transform. Vi fouriertransformerade darfor serien med hjalp av rutinen fft. Absoulutbeloppet av koeficcienterna, vilka ju beskriver frekvensinnehållet, ses i fig XX. Den forsta piken i plotten, $f=535.3$ Hz,  overenstammer ganska bra med den teoretiska grundfrekvensen $c_0/2L\approx 523$ Hz. 
\\Vi testade sedan ett stort hål, s=2, placerat vid x=0.5. Figur XX visar losningen for t=4. Det stora hålet verkar ha samma effekt som att dela flojten på mitten. Vi har en stående våg i varje del.  Vi tittade på losningen i x=0.25 och fouriertransformerade. Frekvensdiagrammet ses i figur XX. Grundfrekvensen verkar nu ha dubblats (spiken ar vid$f=1040$ Hz). Detta ar ju också samma grundfrekvens som om vi valjer $L=0.33/2$, dvs halverar flojtlangden. Vi tittade på olika placeringar, $x_0$, av stora hål och fann att frekvensen for grundsvangningen var ungefar $5523/x_0$.\\
Vi minskade sedan hålstorleken till $s=0.5$. Då hålet placeras vid x=0.5 resp 0.7 fick vi en grundfrekvens på 851.0 resp 729.7. Frekvensdiagrammen visas i fig XX och XX.. Det verkar som att ett mindre hål resulterar i en mindre grundfrekvens.
\begin{figure}[h]
\begin{center}
\scalebox{0.5}{\includegraphics{0205.eps}}
\caption{fft av serien - semilogy}
\end{center}
\end{figure}


\begin{figure}[h]
\begin{center}
\scalebox{0.5}{\includegraphics{stort.eps}}
\caption{Stort hål vid x=0.5}
\end{center}
\end{figure}

\begin{figure}[h]
\begin{center}
\scalebox{0.5}{\includegraphics{0207.eps}}
\caption{fouriertransformen av serien vid x=0.25 - semilogy}
\end{center}
\end{figure}

\subsection{Alternativ berakning av grundfrekvensen}

En stående våg med enveloppen v(x) och frekvensen f kan skrivas $U(t,x)=v(x)exp(i2\pi ft)$(detta var givet i uppgiften), och uppfyller diff-ekvationen \begin{equation}U_{tt}-U_{xx}+d(x)U=0,    U(0)=U(1)=0,\end{equaton} dvs ekvation (8) med G=0. Efter insattning och multiplikation med $exp(-2\pi ft)$ får vi ekvationen for v(x)\begin{equation}-v_{xx}+d(x)v(x)=(2\pi f)^2v(x),    v(0)=v(1)=0,\end{equaton} Vi observerar att ekvation (9) och (10) sager oss att om $u_{tt}-u_{xx}+d(x)u(x)=0,  u(0)=u(1)=0$ så kan vi approximera losningen genom att betrakta ODEn $-\textbf{\emph{u}}_t_t=-(\textbf{\emph{A}}-\textbf{\emph{D}})\textbf{\emph{u}}$. 
\\ Detta ger att $-\textbf{\emph{A}}+\textbf{\emph{D}}$ ar numeriska approximationen av operatorn $\partial_{xx}+d(x)$. Genom att berakna de minsta egenvardena $\lambda_{min}$till $-\textbf{\emph{A}}+\textbf{\emph{D}}$ kan vi alltså uppskatta $f_0$ till $\frac{c_0}{L}\sqrt{\lambda_{min}}/2\pi$ ($L$ fortfarande $0.33$).Vi berknade på detta satt grundfrekvensen, $f0$, till $521.2$ då $s=2$ och $x=0.5$$f_0$ Detta stammer val med vad vi tidigare erhållit. For $s=0.2$ och $x=0.7$ fick vi $f_0=669.6$ och $x=0.5$ gav $f_0=793.3$. Detta ar dock en bit ifrån de $730$ resp $850$ vi fick genom fourieranalysen.



\end{document}