MATLAB 的 bvp5c - 编程之家

如何解决MATLAB 的 bvp5c

我正在尝试解决以下一组二阶非线性和耦合 ODE：

    0 = ?² ℎ″(?) − ? ℎ′(?) + ?² ?(?)² [1−ℎ(?)],0 = ?² ?″(?) + ? ?'(?) − ? ?² ?(?) [?(?)² + ?(?)² − 2],0 = ?² ?″(?) + ? ?′(?) − (1/2) ?(?) [1−ℎ(?)]² − ? ?² ?(?) [?(?)² + ?(?)² − 2].

（很抱歉 EQ 看起来不漂亮。我只习惯于 LaTeX 语法）。好吧，除了方程之外，我还有以下边界条件：

    f'(0) = 0,f'(x → ∞) = 0. I also kNow that f (x → ∞) = 1,h(0)  = 0,h(x → ∞)  = 1,g(0)  = 0,g(x → ∞)  = 1 .

此外，我还希望一阶导数 h'(x)、f'(x) 和 g'(x) 对于某个有限值 x 趋于零，然后保持为零。也就是说，我知道我的解决方案必须达到渐近值并且之后保持不变。换句话说，我知道函数 h(x)、f(x) 和 g(x) 必须“饱和”。

使用 MATLAB，我尝试了以下解决方案：

xmin=1e-3;
xmax=50;
guess = [ 1 1 1 0 0 0];


xmesh = linspace(xmin,xmax,1e5);
solinit = bvpinit(xmesh,guess);%The last vector is my guess.
options = bvpset('RelTol',1e-5,'NMax',5e6); 

sol = bvp5c(@deriv,@bcs,solinit,options);

Sxint = deval(sol,xmesh);
figure
plot(sol.x(1,:),sol.y(1,'k-');
hold on 
plot(sol.x(1,sol.y(2,'m-');
hold on
plot(sol.x(1,sol.y(3,'k-')
hold off
axis([0 xmax -0.2 1.5])

function dydx = deriv(x,y)
lambda=0.5;
dydx= [ y(4) %The vector y() was keeping the following results: y=(h,f,g,h',f',g')
        y(5)
        y(6)
        (1/x)*y(4) - (y(3)^2)*(1-y(1))
        -(1/x)*y(5) + (lambda)*y(2)*(y(2)^2 + y(3)^2 - 2)
        -(1/x)*y(6) + (1/(2*x^2))*y(3)*((1-y(1))^2) + (lambda)*y(3)*(y(2)^2 + y(3)^2 - 2)];
end

% boundary conditions 

function res = bcs(ya,yb)
res = [ya(1)
       yb(1) - 1
       yb(2) - 1
       ya(3)
       yb(3) - 1
       ya(5)];
end

好吧，在复制和粘贴我的代码时我可能会犯一些小错误，这段代码给了我一个解决方案，它只在边界处采用所需的值（“无穷大”）。我可以使用越来越大的 xmax 值，即使如此，解决方案也永远不会在其无穷大值处开始饱和。

我尝试使用更好的猜测，基于这些方程的解析解，x 的值很小，但没有任何东西给我一个好的解决方案。这就是我在这里寻求一些建议的原因。你怎么认为？由于第三个 ODE 中的 1/x²，MATLAB 是否无法解决这个问题？

谢谢！

解决方法

乍一看，对于大 x，您会进入 f,g 一个带有哈密顿量的机械系统

H = 0.5·(f'(x)²+g'(x)²) - 0.25·?·(f(x)²+g(x)²-2)²

摩擦力逐渐下降到零。势能是一座两侧朝向 -oo 的山丘和圆形顶部的山谷，就像一座古老的风化火山。

初始条件表示移动从压痕的中心开始。使用您当前的方法，您可以在有限的时间内到达压痕周围的边缘，因此速度非零。我希望解决方案继续从潜在能源格局的外部下降。

您正在尝试渐近地到达 f=g=1。这对应于能量表面 H=0。直觉（可能是错误的）和必须收敛到零的角动量表明最终方法必须是径向的，即 f(x)=g(x)，因此

f'(x)² = ?·(f(x)²-1)² ==> f'(x) = c·(1-f(x)²),c=sqrt(?)
f(x)=tanh(c*(x-d))

这种渐近行为表明无穷远的边界条件被 x=T 处的条件替换

f'(T)-c*(1-f(T)^2)=0
g'(T)-c*(1-g(T)^2)=0
h'(T)+g(T)*(1-h(T))=0

最后一个杀死 h''+g²·(1-h)=0 的指数增长项。

还需要对 x=0 进行一些合理的近似，不考虑三次项和高阶项，可以得到以下结果：

h(x)=h_2·x² 遵循第一个等式，其中 h_2 免费，
f(x)=f_0+f_2·x² 导致 4*f_2 = ?·f_0·(f_0²-2)，
g(x)=g_1·x+g_2·x² 通向 g_1·x+4·g_2·x² = (g_1·x+g_2·x²)/2，因此 g_1=g_2=0。

接下来的数值实验。或不。到目前为止，上边界的近似值不足以强制收敛。更合理的结果是在角度 0 到 pi/4 的圆段附近