Solving a PDE using spectral methods in Matlab

[nickthequick]

Hi,

I'm trying to solve the PDE
$$u_{tt}-c^2u_{xx}=\omega E_{xx}$$ where

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

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

$\hat{u}_{xx}+\frac{\sigma^2}{c^2} \hat{u}=\frac{\omega}{c^2}\hat{E}_{xx}$

Expanding the RHS as a Fourier series -

$\hat{E}_{xx}=\sum_n a_n \cos k_n x$ means that

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

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 $\hat{E}_{xx}$
**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

 

[NateD]

Unfortunately, there are some limitations to this method that could be contributing to the seemingly incongruous solutions. The most important limitation is that Fourier transforms assume the solution is periodic in both the spatial and temporal domains. If this is not the case, then the Fourier transform will not accurately capture the behavior of the solution and the resulting solution will not be valid. Additionally, this method may also suffer from numerical errors due to discretization of the domain, truncation of the Fourier series, and other numerical approximations. To ensure accuracy, it is important to use a numerical method that considers the boundary conditions and the physical properties of the problem.