I try to solve numerically coupled "sine-Gordon-like" system of the 2nd order differential equations by Matlab
$\frac{{{d^2}\phi }}{{d{\varphi ^2}}} - \left( {1 + \frac{1}{{{\kappa _2}}}} \right)\sin \phi - \sin \theta - \frac{1}{{{\kappa _2}}}\sin \left( {\theta - \phi } \right) = 0,$
$\frac{{{d^2}\theta }}{{d{\varphi ^2}}} - \sin \phi - \left( {1 + \frac{1}{{{\kappa _3}}}} \right)\sin \theta + \frac{1}{{{\kappa _3}}}\sin \left( {\theta - \phi } \right) = 0$,
with the boundary conditions
$\phi \left( {2\pi } \right) - \phi \left( 0 \right) = 2\pi {n_2}$, $\frac{{d\phi }}{{d\varphi }}\left( 0 \right) = \frac{{d\phi }}{{d\varphi }}\left( {2\pi } \right)$,
$\theta \left( {2\pi } \right) - \theta \left( 0 \right) = 2\pi {n_3}$, $\frac{{d\theta }}{{d\varphi }}\left( 0 \right) = \frac{{d\theta }}{{d\varphi }}\left( {2\pi } \right)$,
where $n_2$ and $n_3$ are integer numbers and coefficients $\kappa_2$ and $\kappa_3$ are positive values.
Below is my Matlab code for the solutions of the system
function soliton_3band
tic
options = bvpset('RelTol',10^(-3),'NMax',10^6);
solinit = bvpinit(linspace(0,2*pi,10),@init);
sol = bvp4c(@Eqs,@bc,solinit,options);
xint = linspace(0,2*pi,70000);
Sxint = deval(sol,xint);
plot(xint,Sxint(2,:),'LineWidth',3)
grid on
set(gca,'FontName','Times New Roman','FontSize',24)
set(gca,'xlim',[0 2*pi],'xtick',[0:pi/2:2*pi],'xticklabel',{'0' 'p/2' 'p' '3p/2' '2p'}, 'fontname','symbolpi')
set(gca,'ylim',[-2*pi 2*pi],'ytick',[-2*pi:pi/2:2*pi],...
'yticklabel',{'-2p' '-3p/2' '-p' '-p/2' '0' 'p/2' 'p' '3p/2' '2p'}, 'fontname','symbolpi')
function dydx = Eqs(x,y)
k2=0.024;
k3=2;
dydx = [y(2);
(1+1/k2)*sin(y(1))+sin(y(3))+1/k2*sin(y(3)-y(1));
y(4);
sin(y(1))+(1+1/k3)*sin(y(3))-1/k3*sin(y(3)-y(1))];
end
function res = bc(ya,yb)
n2=1;
n3=0;
res = [ ya(1) - yb(1)+2*pi*n2;
ya(2) - yb(2);
ya(3) - yb(3)+2*pi*n3;
ya(4) - yb(4)];
end
function v = init(x)
v = [ sin(x);
cos(x);
sin(x);
cos(x)];
% v = [ -pi+4*atan(exp(x));
% 4*exp(x)/(1+exp(2*x));
% -pi+4*atan(exp(x));
% 4*exp(x)/(1+exp(2*x))];
% v = [-pi+4*atan(exp(x)); 4*exp(x)/(1+exp(2*x));sin(x);cos(x)];
end
toc
end
For some limiting cases I can get analytical solutions of that system, for instance, for $\kappa_2$ = 0 the system transforms to one algebraic equation and one so-called double sine-Gordon equation for the function $\theta$, which solution in turn can be expressed via Jacobi elliptic functions.
I'd like to compare my analytical solutions for unknown functions with the numerical one with very small (vanishing) $\kappa_2$.
And here I have a problem. I solved the system for some parameters with the decreasing value of $\kappa_2$ (see figure) for the given set of other parameters.
But for some small coefficients $\kappa_2$ Matlab is not able to find solution for $\phi$ because unable to meet the tolerance (even with a larger mesh grid). How can I change my Matlab code to solve the system numerically for small $\kappa$ (for instance for 0.001)?
Btw Comsol can't solve that system also for small $\kappa_2$.
I will appreciate for any help.