- 53

- 0

I'm trying to solve the PDE

[tex] u_{tt}-c^2u_{xx}=\omega E_{xx} [/tex] where

[tex]E= A(x)Sech^2\left(\frac{x-ct}{B(x)}+\delta(x)\right) [/tex]

To this end, I take the temporal Fourier Transform of the PDE to find that the governing equation is the forced simple harmonic oscillator:

[itex] \hat{u}_{xx}+\frac{\sigma^2}{c^2} \hat{u}=\frac{\omega}{c^2}\hat{E}_{xx} [/itex]

Expanding the RHS as a Fourier series -

[itex] \hat{E}_{xx}=\sum_n a_n \cos k_n x [/itex] means that

[tex] \hat{u}=\frac{\omega}{c^2}\sum_n \frac{a_n}{(\sigma/c)^2 -k_n^2}[/tex] and finally one can find u by taking the inverse Fourier transform of [itex]\hat{u}[/itex].

My issue is coding this up, since it appears that I'm getting nothing but nonsense for a solution. I am definitely a novice in numerical methods and I would like some input from someone with some experience about whether or not this type of method is good in this situation. Also, since I believe my code is doing what I intend it to, are there any limitations to this method that could account for my seemingly incongruous solutions?

Thanks!

Nick

**here F3 is the prescribed forcing function and E_hat_t is [itex] \hat{E}_{xx}[/itex]

**positiveFFT finds the positive FT as well as the associated frequencies

% E_hat_t=zeros(length(x1),length(t)/2);

% f_E_t=zeros(length(x1), length(t)/2);

% for i = 1:length(x1);

% [E_hat_t(i,:),f_E_t(i,:)]=positiveFFT(F3(i,1:length(t)),Xs1);

% end

%

% E_hat_x1= zeros((length(x1))/2,length(t)/2);

% f_E_x1=zeros((length(x1))/2,length(t)/2);

% for i =1:(length(EDa(1,)/2

% [E_hat_x1(:,i),f_E_x1(:,i)]=positiveFFT(E_hat_t(:,i),Xs1);

% end

%

u_hat=zeros((length(x1))/2,length(t)/2);

for i=2:(length(x1))/2

for j=2:length(t)/2

u_hat(i,j)=real(omega_c/(g*H).*(E_hat_x1(i,1)+2*sum(E_hat_x1(i,:).*exp(sqrt(-1)*2*pi*f_E_x1(i,:)*x1(i))./((2*pi*f_E_t(i,j)).^2/(g*H)-(2*pi*f_E_x1(i,.^2))));

end

end

% u=zeros(length(x1)/2,length(t)/2);

% for i=2:(length(x1)-1)/2

% for j=2:length(t)/2

% u(i,j)= u_hat(i,1)+2*sum(u_hat(i,:).*exp(sqrt(-1)*f_E_t(i,:)*t(j)));

% end

% end