Weakly nonlinear theories in electrohydrodynamics

[hunt_mat]

I derived an equation describing the free surface of an electrified fluid. I am currently seeking traveling wave solutions for this problem, the equation I am looking at is (1-F) f+\frac{1}{90}h^{4}f^{(4)}+\frac{3}{4h}f^{2}-\frac{1}{2}\Bigg( B-\frac{1}{3}\Bigg) hf''+\frac{Ah}{2}\mathcal{H}(f')+\frac{p}{2\rho g}=0

The \mathcal{H} term is the Hilbert transform, F,h,B,\rho ,g are constants

I have written a numerical program which solves the above equation in MATLAB and it seems to give the sort of solution I am looking for. In order to text this result, if I ignore the nonlinear term, I can derive an analytical solution by the use of Fourier transforms, so if I write:

f=\int_{\mathbb{R}}\alpha (k)e^{ikx}dk,\quad p=\int_{\mathbb{R}}\beta (k)e^{ikx}dk.

I can get the solution

\alpha (k)=-\frac{\beta (k)}{2\rho g}\bigg[ 1-\frac{c}{c_{0}}+\frac{h^{4}k^{4}}{90}+\frac{1}{2}\Bigg( B-\frac{1}{3}\Bigg) h^{2}k^{2}-\frac{Ah|k|}{2}\Bigg]^{-1}

If I compare the above solution to my numerical solution with the nonlinear term ignored, they should, in theory, get the same answer. They don't however. Can anyone see why this should be? The MATLAB program I use is:

L=201; %Defines the number of steps, must always be an integer
x=linspace(-6,6,L); %Defines x
p=exp(-x.^2); %The pressure distribution
f=zeros(1,L); %Initial guess for f is f=0;
Q=0.1; %This the the ratio c/c_0
%The way forward is to use Newton's method to reduce the equation down to a linear algebra problem and then use Octave
%to solve the linear algrbra problem.

%The first step is to compute the Jacobian matrix
%Do this by diferencing between to close points and dividing through by that difference
df=10^-10;
J=zeros(L,L); %This will be the Jacobian.
it=3; %This will be the number of iterations required to compute the solution
X=zeros(L,1);
for k=1:it
for i=1:L
for j=1:L
T=zeros(1,L); %Set all the entries to zero again
T(i)=df;
J(i,j)=(g(x,f+T,p,Q,j)-g(x,f,p,Q,j))/df; %Calculates the (i,j)th element of the Jacobian
end
end

for m=1:L
X(m)=-g(x,f,p,Q,m);
end

%Now solve the system of equations
sol=J\X;
%sol=inv(J)*X;
f=f+sol';

end

%Now plot the solution!
plot(x,f);
xlabel('x');
ylabel('f(x-ct)');

The function g is:

function r=g(x,f,p,Q,i)
h=1;
B=0.2;
E_b=0.1;
g=9.8065;
rho=1;
a_1=1-Q;
a_2=h^4/90;
a_3=0.75/h;
a_4=-0.5*(B-1/3)*h^2;
a_5=0.5*E_b*h;
a_6=1/(2*g*rho);
dx=x(2)-x(1);
f_ghost=[0 0 f 0 0]; %Add in the ghost cells for f
N=length(x); %Prepare for the Hilbert transform
x_half=zeros(1,N-1);
f_half=zeros(1,N-1);
for k=1:N-1
x_half(k)=0.5*(x(k)+x(k+1));
f_half(k)=0.5*(f(k)+f(k+1));
end
grad=gradient(f_half,dx);
hilbert=trapz(x_half,grad./(x_half-x(i)));
%hilbert=0; %Start the Hilbert transform at zero
%for k=2:N-1
%hilbert=hilbert+(f(k+1)-f(k))/(0.5*(x(k+1)+x(k))-x(i));
%end
%Compute the Hilbert transform

%Now comes the formulae for the finite differencing part of the algorithm.
b_1=a_1*f_ghost(i+2);
b_2=(a_2/dx^4)*(f_ghost(i+4)-4*f_ghost(i+3)+6*f_ghost(i+2)-4*f_ghost(i+1)+f_ghost(i)); %Fourth order derivative
b_3=0; %a_3*f_ghost(i+2)^2; %nonlinear term
b_4=(a_4/dx^2)*(f_ghost(i+3)-2*f_ghost(i+2)+f_ghost(i+1)); %Second order derivative
b_5=a_5*hilbert; %Hilbert transform term
b_6=a_6*p(i); %Pressure term
r=b_1+b_2+b_3+b_4+b_5+b_6;

Any idea anyone?
Mat