Блог

Я пытаюсь решить следующий набор нелинейных и связанных ОДУ 2-го порядка:

     0 = 𝑥² ℎ″(𝑥) − 𝑥 ℎ′(𝑥)   𝑥² 𝑔(𝑥)² [1−ℎ(𝑥)], 
    0 = 𝑥² 𝑓″(𝑥)   𝑥 𝑓'(𝑥) − 𝜆 𝑥² 𝑓(𝑥) [𝑓(𝑥)²   𝑔(𝑥)² − 2], 
    0 = 𝑥² 𝑔″(𝑥)   𝑥 𝑔′(𝑥) − (1/2) 𝑔(𝑥) [1−ℎ(𝑥)]² − 𝜆 𝑥² 𝑔(𝑥) [𝑓(𝑥)²   𝑔(𝑥)² − 2].
 

(Мне жаль, что эквалайзеры не выглядят красиво. Я привык только к синтаксису 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
 

Ну, вплоть до некоторых незначительных опечаток, которые я мог бы сделать при копировании и вставке моего кода, этот код дает мне решение, которое принимает только желаемое значение («бесконечность») на границе. Я могу использовать все большие и большие значения my xmax , и даже в этом случае решение никогда не начинает насыщаться по своему значению на бесконечность.

Я попытался использовать лучшее предположение, основанное на аналитическом решении этих уравнений для небольшого значения x , но ничто не дает мне хорошего решения. И именно по этой причине я прошу здесь совета. Что вы думаете? Является ли MATLAB неспособным решить это из 1/x² -за третьего ODE?

Спасибо!

На первый взгляд, для больших 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