Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Solving a PDE using spectral methods in Matlab

  1. Nov 20, 2011 #1
    Hi,

    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
     
  2. jcsd
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Can you offer guidance or do you also need help?
Draft saved Draft deleted