szw3.m

来自「如何利用bvp4c解决微分方程的边值问题的几个例子。」· M 代码 · 共 65 行

function ex1bvp
%The problem is
%   
%      v'' =  a*v-v*w
%      w'' = (1/b^2)*(w-0.5*v^2)
%
%   The interval is [-15 15] and the boundary conditions are
%   
%      w(-15) = v(-15) = 0,  w(15) = v(15) = 0,
global a b
a=1;
b=1;
solinit = bvpinit(linspace(-20,20,11),@ex1init);
options = bvpset('Stats','on','RelTol',1e-9);

sol = bvp4c(@ex1ode,@ex1bc,solinit,options,a,b);

% The solution at the mesh points
x = sol.x;
y = sol.y;

clf reset
plot(x,y(1,:),'b-')
axis auto
xlabel('x')
shg
hold on
%clf reset
plot(x,y(3,:),'r-')
axis auto
xlabel('x')
shg
% --------------------------------------------------------------------------

function dydx = ex1ode(x,y,a,b)

%   The components of y correspond to the original variables
%   as  y(1) = v, y(2) = v', y(3) = w, y(4) = w'.

dydx =  [ y(2)
         a*y(1)-y(3)*y(1)
         y(4)
         1/b^2*(y(3)-0.5*y(1)^2) ];

%-------------------------------------------------------------------------

function res = ex1bc(ya,yb,a,b)

res = [ ya(1)
        yb(1)
        ya(2)
        yb(2)];

%-------------------------------------------------------------------------

function v = ex1init(x)
%EX1INIT guess for problem of the BVP4C.
global a

v = [2*sqrt(a)*sech(sqrt(a)*x) 
     2*sqrt(a)*sech(sqrt(a)*x)*tanh(sqrt(a)*x)
     2*a*(sech(sqrt(a)*x))^2
     -4*a*(sech(sqrt(a)*x))^2*tanh(sqrt(a)*x)];