I tried to run it on another computer, previously it ran and converged within 4 iteration to an error of 10^-8, now it won't converge for some reason.
The code:
clear;
L=201; %Defines the number of steps, must always be an integer
x=linspace(-6,6,L); %Defines x
B=0.2; %This is the Bond number.
p=0.25*B*exp(-x.^2); %The pressure distribution
E_b=0.2; %This is the electric Bond number.
f=zeros(1,L+1);
f(1,1:L)=0.01*exp(-x.^2); %have some small initial guess for the wave profile
choice=false; %This selects the option to fix the Froude number of the minimum of the wave.
Q=-0.6; %This is the Froude number (choice=true), or the minimum of the profile (choice=false)
f(L+1)=Q; %this is the Foude number, either an exact value or an initial guess.
%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+1,L+1); %This will be the Jacobian.
X=zeros(L+1,1);
error=0.5;
while (error>10^-10)
for i=1:L+1
for j=1:L+1
T=zeros(1,L+1); %Set all the entries to zero again
T(j)=df;
J(i,j)=(g(x,f+T,p,Q,i,choice,B,E_b)-g(x,f-T,p,Q,i,choice,B,E_b))/(2*df); %Calculates the (i,j)th element of the Jacobian
end
end
for m=1:L+1
X(m)=-g(x,f,p,Q,m,choice,B,E_b);
end
%Now solve the system of equations
sol=J\X;
error=norm(X)
f=f+sol';
end
%Now plot the solution!
plot(x,f(1:L));
xlabel('x');
ylabel('f(x-ct)');
Using:
function r=g(x,f,p,Q,i,flag,B,E_b)
N=length(x);
h=1;
a_1=1-f(N+1);
a_2=h^4/90;
a_3=1.5/h;
a_4=-0.5*(B-1/3)*h^2;
a_5=0.5*E_b*h;
dx=x(2)-x(1);
ff=f(1:N);
f_ghost=[0 0 ff 0 0]; %Add in the ghost cells for f
%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*(ff(k)+ff(k+1));
end
grad=gradient(f_half,dx);
if (i<N+1)
hilbert=trapz(x_half,grad./(x_half-x(i)));
end
%Now comes the formulae for the finite differencing part of the algorithm.
if (i<N+1)
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=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=p(i); %Pressure term
r=b_1+b_2+b_3+b_4+b_5+b_6;
else
if (flag==true)
r=f(N+1)-Q;
else
r=f(1+0.5*(N-1))-Q;
end
end
