Нелинейное дифференциальное уравнение, как мне решить это численно в MATLAB?

Question

Я работал над проектом, в котором мне нужно найти решение данного нелинейного дифференциального уравнения, см. Рисунок ниже:

Теперь я попытался использовать matlabsвстроенная функция bvp4c, но синтаксис сложен, и я не знаю, надежны ли результаты.Для некоторых значений функция bvp4c генерирует только ошибку.Я также должен принять во внимание граничные условия, см. Рисунки ниже:

Прошу прощения за ужасные размеры фигур.Теперь я знаю, что это не математический форум, но мне нужно решить это численно, и я хочу решить это наилучшим доступным методом.Мой код на данный момент показан ниже:

function [theta_0 y2]=flow_BVP

theta_0=linspace(0,1,1000); % pi/18
solinit2=bvpinit(theta_0,[0 1 1]);
sol2=bvp4c(@flow_ode,@flow_bc,solinit2);
x2=sol2.x;
y2=sol2.y(1,:);


hold on
%plot(x1,y1) %gammal
plot(x2,y2) %ny
%hold off


function v=flow_init(x)
v=[sin(x); 1; 1];

function dydx=flow_ode(x,y)
q=0.0005;
v=1;
dydx = [y(2); y(3); 2*q/v*y(1)*y(2)-4*y(2)];

function res=flow_bc(ya,yb)
res=[ya(1);yb(1);ya(2)-5.59];

Чтобы повторить мой вопрос, какой метод самый лучший, самый простой, простой для понимания и реализации для решения этой проблемы?Возможно, стрельба?

С наилучшими пожеланиями SimpleP.

Редактировать Результат, который я получил до сих пор, это

Сюжет показывает f против \ theta.Интеграл равен 1, как и должно быть.

Answer 1

Общий способ включения интеграла состоит в добавлении антипроизводного F из f в систему ODE.То есть, поскольку четвертый компонент и переменная

F' = f  with   F(0)=0,  F(alpha)=1

и другие компоненты должны быть сдвинуты на один индекс,

function v=flow_init(x)
v=[sin(x); 1; 1; 1-cos(x)];

function dydx=flow_ode(x,y)
% y is [ f, f', f'', F ]
q=0.0005;
v=1;
dydx = [y(2); y(3); 2*q/v*y(1)*y(2)-4*y(2); y(1)];

function res=flow_bc(ya,yb)
res=[ya(1); ya(4); yb(1); yb(4)-1];

---

Использование python:

q, v = 0.0005, 1
def flow_ode(t,u): return [ u[1], u[2], 2*q/v*u[0]*u[1]-4*u[1], u[0] ]

def flow_bc(u0, u1): return [ u0[0], u0[3], u1[0], u1[3]-1 ]

x = x_init = np.linspace(0,1,11);
u_init = [ 6*x*(1-x), 0*x, 0*x, x ]

res = solve_bvp(flow_ode, flow_bc, x_init, u_init, tol = 1e-5)

print res.message
if res.success:
     x = x_sol = np.linspace(0,1,201);
     u_sol = res.sol(x_sol);
     plt.subplot(2,1,1)
     plt.plot(x_sol, u_sol[0]); plt.plot(x, 6*x*(1-x), lw=0.5); plt.grid()
     plt.subplot(2,1,2)
     plt.plot(x_sol, u_sol[3]); plt.grid()
     plt.show()

Как видно, первоначальное предположение здесь довольно близко.Поскольку ОДУ представляет собой небольшое возмущение, равное 4f'+f'''=0, решение должно быть близко к его решению a+b*sin(2x)+c*cos(2x), которое с граничными условиями оценивается как

f(x)=A * [ (1-cos(2))*sin(2*x)-sin(2)*(1-cos(2*x)) ]
    = 4*A*sin(1) * sin(x)*sin(1-x)

с A, так что интеграл равен единице..

---

Если бы значения параметров в ODE были переключены, q=1 и v=0.0005, решение имело бы граничные слои размером sqrt(v/q).