Решение краевой задачи (bvp4c matlab)

@Rose12345 · Регистрация: 09.07.2016

Author24 — интернет-сервис помощи студентам

Имеется система дифференциальных уравнений. Необходимо решить краевую задачу - заданы некоторые начальные и некоторые конечные значения. Использую функцию bvp4c, задача решается, получаю набор начальных значений. НО! Если подставить полученные начальные значения и проинтегрировать систему ДУ с помощью ode45, результат получается другим. Код программы прилагаю. В чем может быть ошибка?

Кликните здесь для просмотра всего текста

Matlab M

clear, clc, cla
 
global mu_ P_ We_ m0_ Rz_ Vz_ Rm_ Vm_
mu_ = 1;
P_ = 0.03;
We_ = 0.5268;
m0_ = 1;
 
Rz_ = [1.002972 0.018841 -0.7333937*(10^-5)];
Vz_ = [-0.035064 0.996249 -0.653687*(10^-4)];
 
Rm_ = [1.183485 -0.712992 -0.0439577];
Vm_ = [0.4505997 0.766424 0.00501665];
 
N = 10801;
Tp_ = 7.1675;
tin = linspace(0, Tp_, N);
yin = [Rz_(1) Rz_(2) Rz_(3) Vz_(1) Vz_(2) Vz_(3) 10^-10 10^-10 10^-10 10^-10 10^-10 10^-10 m0_];
init = bvpinit(tin, yin);
sol = bvp5c(@odefunc, @bcfunc, init);
x  = sol.x;
y1 = sol.y(1,:);
y2 = sol.y(2,:);
y3 = sol.y(3,:);
y4 = sol.y(4,:);
y5 = sol.y(5,:);
y6 = sol.y(6,:);
y7 = sol.y(7,:);
y8 = sol.y(8,:);
y9 = sol.y(9,:);
y10 = sol.y(10,:);
y11 = sol.y(11,:);
y12 = sol.y(12,:);
y13 = sol.y(13,:);
 
%% Начальные и конечные значения
Yab = vpa([sol.y(:,1) sol.y(:,end)],10)
 
%% Проверка
Y0 = [Rz_(1) Rz_(2) Rz_(3) Vz_(1) Vz_(2) Vz_(3) y7(1) y8(1) y9(1) y10(1) y11(1) y12(1) m0_];
t = 0:(6.6360e-04):Tp_;
[T,Y] = ode23(@(t,y)odefunc(t,y),t,Y0);

Matlab M

function dy = odefunc(t,y)
global mu_ P_ We_ %m0_ Rz_ Vz_ Rm_ Vm_
 
dy = zeros(13,1);
dy(1) = y(4);
dy(2) = y(5);
dy(3) = y(6);
dy(4) = -mu_*y(1)*((y(1)^2+y(2)^2+y(3)^2)^(-3/2))+P_/y(13)*y(7)*((y(7)^2+y(8)^2+y(9)^2)^(-1/2));
dy(5) = -mu_*y(2)*((y(1)^2+y(2)^2+y(3)^2)^(-3/2))+P_/y(13)*y(8)*((y(7)^2+y(8)^2+y(9)^2)^(-1/2));
dy(6) = -mu_*y(3)*((y(1)^2+y(2)^2+y(3)^2)^(-3/2))+P_/y(13)*y(9)*((y(7)^2+y(8)^2+y(9)^2)^(-1/2));
dy(7) = -y(10);
dy(8) = -y(11);
dy(9) = -y(12);
dy(10) = y(7)*mu_*((y(1)^2+y(2)^2+y(3)^2)^(-3/2))-3*mu_*y(1)*((y(1)^2+y(2)^2+y(3)^2)^(-5/2))*(y(7)*y(1)+y(8)*y(2)+y(9)*y(3));
dy(11) = y(8)*mu_*((y(1)^2+y(2)^2+y(3)^2)^(-3/2))-3*mu_*y(2)*((y(1)^2+y(2)^2+y(3)^2)^(-5/2))*(y(7)*y(1)+y(8)*y(2)+y(9)*y(3));
dy(12) = y(9)*mu_*((y(1)^2+y(2)^2+y(3)^2)^(-3/2))-3*mu_*y(3)*((y(1)^2+y(2)^2+y(3)^2)^(-5/2))*(y(7)*y(1)+y(8)*y(2)+y(9)*y(3));
dy(13) = -P_/We_;
 
end

Matlab M

function dy = bcfunc(y0,yN)
global  m0_ Rz_ Vz_ Rm_ Vm_ %mu_ P_ We_
dy = zeros(5,1);
dy(1) = y0(1)-Rz_(1);
dy(2) = y0(2) - Rz_(2);
dy(3) = y0(3) - Rz_(3);
 
dy(4) = y0(4)-Vz_(1);
dy(5) = y0(5) - Vz_(2);
dy(6) = y0(6) - Vz_(3);
 
dy(7) = yN(1)-Rm_(1);
dy(8) = yN(2) - Rm_(2);
dy(9) = yN(3) - Rm_(3);
 
dy(10) = yN(4)-Vm_(1);
dy(11) = yN(5) - Vm_(2);
dy(12) = yN(6) - Vm_(3);
dy(13) = y0(13)-m0_;
 
end

Так выглядит решение при использовании bvp4c и ode45

@Centurio · 28.06.2020, 16:52

Rose12345, ode45 не используется для решения краевых задач.

@Rose12345 · 28.06.2020, 17:07 **[ТС]**

Centurio, С помощью ode45 я проинтегрировала уравнения, используя начальные значения, которые получила из bvp4c

@Centurio · 29.06.2020, 11:36

Rose12345, вы написали бы лучше саму систему уравнений и краевые условия. А то по коду разбирать, что вы решаете, тяжело, а если код у вас ещё и с ошибкой, то вообще бессмысленно.

@Rose12345 · 03.07.2020, 15:46 **[ТС]**

Имеется система из 13 ДУ.
P = 0.03; W = 0.53. Время интегрирования t = 7.17.
В начальный момент времени известны: y1(0) = 1.002972; y2(0) = 0.018841; y3(0) = -0.7333937*(10^-5);
y4(0) = -0.035064; y5(0) =0.996249; y6(0) = -0.653687*(10^-4); y13(0) = 1. Необходимо найти такие y7(0), y8(0) ,y9(0), y10(0), y11(0), y12(0), чтобы выполнялись конечные условия:

y1(t) = 1.183485; y2(t) = -0.712992; y3(t) = -0.0439577;
y4(t) = 0.4505997; y5(t) =0.766424; y6(t) = 0.00501665.
Как это можно сделать? Какие функции использовать?

$https://www.cyberforum.ru/cgi-bin/latex.cgi?\dot{y_1}=y_4; \ \\\dot{y_2}=y_5;\\\dot{y_3}=y_6;\\\dot{y_4}=-\frac{y_1}{(sqrt(y_1^2+y_2^2+y_3^2))^3}+\frac{P}{y_1_3}*\frac{y_7}{sqrt(y_7^2+y_8^2+y_9^2)};\\\dot{y_5}=-\frac{y_2}{(sqrt(y_1^2+y_2^2+y_3^2))^3}+\frac{P}{y_1_3}*\frac{y_8}{sqrt(y_7^2+y_8^2+y_9^2)};\\\dot{y_6}=-\frac{y_3}{(sqrt(y_1^2+y_2^2+y_3^2))^3}+\frac{P}{y_1_3}*\frac{y_9}{sqrt(y_7^2+y_8^2+y_9^2)};\\\dot{y_7} = -y_1_0;\\\dot{y_8} = -y_1_1;\\\dot{y_9} = -y_1_2;\\\dot{y_1_0} = \frac{y_7}{(sqrt(y_1^2+y_2^2+y_3^2))^3}-\frac{3*y_1*(y_7*y_1+y_8*y_2+y_9*y_3)}{(sqrt(y_1^2+y_2^2+y_3^2))^5};\\\dot{y_1_1} = \frac{y_8}{(sqrt(y_1^2+y_2^2+y_3^2))^3}-\frac{3*y_2*(y_7*y_1+y_8*y_2+y_9*y_3)}{(sqrt(y_1^2+y_2^2+y_3^2))^5};\\\dot{y_1_2} = \frac{y_9}{(sqrt(y_1^2+y_2^2+y_3^2))^3}-\frac{3*y_3*(y_7*y_1+y_8*y_2+y_9*y_3)}{(sqrt(y_1^2+y_2^2+y_3^2))^5};\\\dot{y_1_3} = -\frac{P}{W}$

@nuHrBuH · 05.07.2020, 13:25

Нужно использовать bvp4c как сказал выше Centurio
Простой пример:
Уравнение: $https://www.cyberforum.ru/cgi-bin/latex.cgi?m\ddot x -c \dot x + k x = 0$
Перепишем его в виде:
$https://www.cyberforum.ru/cgi-bin/latex.cgi?\dot x_1 = x_2$
$https://www.cyberforum.ru/cgi-bin/latex.cgi?\dot x_2 = (c x_2-k x_1)/m$
Граничные и начальтные условия: $https://www.cyberforum.ru/cgi-bin/latex.cgi?y(0)=0, \dot y(end)=1$

Matlab M

c = 0.8;
k = 1.4;
m = 1;
fode = @(t,y) [x(2); (c*x(2)-k*x(1))/m];    % дифф уравнение
fbound = @(xa,xb) [xa(1);xb(2)-1];          % граничные/начальные условия
 
init = bvpinit([0 5],[1 0]);                % начальное предположение
s = bvp4c(fode,fbound,init);                % решай
 
plot(s.x,s.y)                               % рисуй

@Rose12345 · 05.07.2020, 14:24 **[ТС]**

nuHrBuH, если использовать тот код, который вы прислали, то используя bvp4c получили решение краевой задачи (графики) и начальные условия (s.y(:,1) = [0 0.2607]) . Так как это система ДУ, то зная начальные условия (s.y(:,1) = [0 0.2607]), можно проинтегрировать систему ДУ и тоже получить решение. Если можно так сказать, выполнить проверку, правильно ли использовали bvp4c. Если воспользоваться полученным приближением - (s.y(:,1) = [0 0.2607]) и использовать ode45 для интегрирования ДУ - fode на интервале [0 5] как и прежде, то в Вашем случае получится, что графики "сойдутся" - значит, что bvp4c использовали верно и верно нашли начальное приближение и, следовательно, нашли верное решение краевой задачи (прикрепляю полученные графики). В моем коде (самое первое сообщение) я решила систему ДУ (краевую задачу), используя bvp4c, нашла начальные значения и, используя эти значения как начальные условия для интегрирования уравнений (с помощью ode45), получила совершенно другое решение. В моем случае, графики не сошлись. Так как вероятность совершить ошибку в использовании ode 45 гораздо ниже, чем в bvp4c, то скорее всего ошиблась в bvp4c. Я делала по аналогии с Вашим примером. Но результат другой. Мне как раз нужно получить начальные значения, т.е. y7,y8,y9,y,10,y11,y12. Где в моем коде может быть ошибка? Или может быть в матлабе есть еще функции, позволяющие найти "неизвестные" начальные условия?

@nuHrBuH · 05.07.2020, 19:54

ВОт мои рекомендации:
рекомендация 1

Matlab M

1 2	% yin = [Rz_(1) Rz_(2) Rz_(3) Vz_(1) Vz_(2) Vz_(3) 10^-10 10^-10 10^-10 10^-10 10^-10 10^-10 m0_]; yin = [Rz_ Vz_ 1e-10*ones(1,6) m0_];

рекомендация 2

Matlab M

function dy = odefunc(t,y)
global mu_ P_ We_ %m0_ Rz_ Vz_ Rm_ Vm_
 
dy = zeros(13,1);
dy(1:3) = y(4:6);
y123 = sqrt(sum(y(1:3).^2));
y789 = sqrt(sum(y(7:9).^2));
dy(4:6) = -mu_*y(1:3)/y123^3 + P_/y(13)*y(7:9)/y789;
dy(7:9) = -y(10:12);
C1 = sum(y(7:9).*y(1:3));
dy(10:12) = y(7:9)*mu_/y123^3 - 3*mu_*y(1:3)*C1/y123^5;
dy(13) = -P_/We_;
 
end

Действительно решения не совпадают
Пробовал уменьшать шаг интеграции, но безуспешно. Может кто-нибудь из экспертов заинтересуется?

@Centurio Модератор 1700 / 1552 / 520 Регистрация: 13.09.2015 Сообщений: 5,370
	28.06.2020, 16:52	2
	Rose12345, ode45 не используется для решения краевых задач. 0

@Rose12345 0 / 0 / 0 Регистрация: 09.07.2016 Сообщений: 7
	28.06.2020, 17:07 [ТС]	3
	Centurio, С помощью ode45 я проинтегрировала уравнения, используя начальные значения, которые получила из bvp4c 0

@Centurio Модератор 1700 / 1552 / 520 Регистрация: 13.09.2015 Сообщений: 5,370
	29.06.2020, 11:36	4
	Rose12345, вы написали бы лучше саму систему уравнений и краевые условия. А то по коду разбирать, что вы решаете, тяжело, а если код у вас ещё и с ошибкой, то вообще бессмысленно. 0

@Rose12345 0 / 0 / 0 Регистрация: 09.07.2016 Сообщений: 7
	03.07.2020, 15:46 [ТС]	5
	Имеется система из 13 ДУ. P = 0.03; W = 0.53. Время интегрирования t = 7.17. В начальный момент времени известны: y1(0) = 1.002972; y2(0) = 0.018841; y3(0) = -0.7333937(10^-5); y4(0) = -0.035064; y5(0) =0.996249; y6(0) = -0.653687(10^-4); y13(0) = 1. Необходимо найти такие y7(0), y8(0) ,y9(0), y10(0), y11(0), y12(0), чтобы выполнялись конечные условия: y1(t) = 1.183485; y2(t) = -0.712992; y3(t) = -0.0439577; y4(t) = 0.4505997; y5(t) =0.766424; y6(t) = 0.00501665. Как это можно сделать? Какие функции использовать? $https://www.cyberforum.ru/cgi-bin/latex.cgi?\dot{y_1}=y_4; \ \\\dot{y_2}=y_5;\\\dot{y_3}=y_6;\\\dot{y_4}=-\frac{y_1}{(sqrt(y_1^2+y_2^2+y_3^2))^3}+\frac{P}{y_1_3}\frac{y_7}{sqrt(y_7^2+y_8^2+y_9^2)};\\\dot{y_5}=-\frac{y_2}{(sqrt(y_1^2+y_2^2+y_3^2))^3}+\frac{P}{y_1_3}\frac{y_8}{sqrt(y_7^2+y_8^2+y_9^2)};\\\dot{y_6}=-\frac{y_3}{(sqrt(y_1^2+y_2^2+y_3^2))^3}+\frac{P}{y_1_3}\frac{y_9}{sqrt(y_7^2+y_8^2+y_9^2)};\\\dot{y_7} = -y_1_0;\\\dot{y_8} = -y_1_1;\\\dot{y_9} = -y_1_2;\\\dot{y_1_0} = \frac{y_7}{(sqrt(y_1^2+y_2^2+y_3^2))^3}-\frac{3y_1(y_7y_1+y_8y_2+y_9y_3)}{(sqrt(y_1^2+y_2^2+y_3^2))^5};\\\dot{y_1_1} = \frac{y_8}{(sqrt(y_1^2+y_2^2+y_3^2))^3}-\frac{3y_2(y_7y_1+y_8y_2+y_9y_3)}{(sqrt(y_1^2+y_2^2+y_3^2))^5};\\\dot{y_1_2} = \frac{y_9}{(sqrt(y_1^2+y_2^2+y_3^2))^3}-\frac{3y_3(y_7y_1+y_8y_2+y_9y_3)}{(sqrt(y_1^2+y_2^2+y_3^2))^5};\\\dot{y_1_3} = -\frac{P}{W}$ 0

@Rose12345 0 / 0 / 0 Регистрация: 09.07.2016 Сообщений: 7
	05.07.2020, 14:24 [ТС]	7
	nuHrBuH, если использовать тот код, который вы прислали, то используя bvp4c получили решение краевой задачи (графики) и начальные условия (s.y(:,1) = [0 0.2607]) . Так как это система ДУ, то зная начальные условия (s.y(:,1) = [0 0.2607]), можно проинтегрировать систему ДУ и тоже получить решение. Если можно так сказать, выполнить проверку, правильно ли использовали bvp4c. Если воспользоваться полученным приближением - (s.y(:,1) = [0 0.2607]) и использовать ode45 для интегрирования ДУ - fode на интервале [0 5] как и прежде, то в Вашем случае получится, что графики "сойдутся" - значит, что bvp4c использовали верно и верно нашли начальное приближение и, следовательно, нашли верное решение краевой задачи (прикрепляю полученные графики). В моем коде (самое первое сообщение) я решила систему ДУ (краевую задачу), используя bvp4c, нашла начальные значения и, используя эти значения как начальные условия для интегрирования уравнений (с помощью ode45), получила совершенно другое решение. В моем случае, графики не сошлись. Так как вероятность совершить ошибку в использовании ode 45 гораздо ниже, чем в bvp4c, то скорее всего ошиблась в bvp4c. Я делала по аналогии с Вашим примером. Но результат другой. Мне как раз нужно получить начальные значения, т.е. y7,y8,y9,y,10,y11,y12. Где в моем коде может быть ошибка? Или может быть в матлабе есть еще функции, позволяющие найти "неизвестные" начальные условия? Миниатюры 0

	05.07.2020, 19:54