Как реализовать шесть начальных условий для системы двух связанных дифференциальных уравнений второго порядка? - PullRequest
Новая
       81

Как реализовать шесть начальных условий для системы двух связанных дифференциальных уравнений второго порядка?

0 голосов
/ 24 октября 2018

Я новичок в кодировании.Основная информация о моей проблеме: r1 и r2 - две переменные;u1 = dr1 / dt, u2 = dr2 / dt и du1 / dt = d ^ 2r1 / dt ^ 2, du2 / dt = d ^ 2r2 / dt ^ 2.В коде Matlab: r (1) подразумевает r1, r (2) -> u1, r (3) -> r2, r (4) -> u2.rdot (2) - это выражение для du1 / dt, а rdot (4) - выражение для du2 / dt.В идеале мне нужно всего 4 начальных условия: r1 (0), u1 (0), r2 (0), u2 (0), которые равны 10d-6, 0, 5d-6, 0. Но в моем случае du1 / dtимеет зависимость от du2 / dt и наоборот.Смотрите последний срок T1_1 и T2_1.И идеальная микросхема для du1 / dt и du2 / dt равна 0. Но как мне реализовать это в моем коде?

Мой код здесь. 

function rdot = f(t, r)

P_stat = 1.01325d5;
P_v = 2.3388d3;
mu = 1.002d-3;
sigma = 72.8d-3;
c_s = 1481d0;
poly_exp = 1.4d0;
rho = 998.2071d0;
f_s = 20d3;
P_s = 1.01325d5;
r1_eq = 10d-6;
r2_eq = 4d-6;
d = 1d-3;


rdot(1) = r(2);
P1_bw = ( (P_stat - P_v + (2.d0*sigma/r1_eq))*((r1_eq/r(1))^(3.d0*poly_exp)) ) - (2.d0*sigma/r(1)) - (4.d0*mu*r(2)/r(1));
P1_ext = P_s*sin(2.d0*pi*f_s*(t + (r(1)/c_s)));
T2_1 = ((2.d0*r(3)*(r(4)^2.d0)) + ((r(3)^2.d0)*rdot(4)))/d;
T2_4 = (1.d0 - (r(2)/c_s))*r(1);
T2_5 = 1.5d0*(1.d0 - (r(2)/(3.d0*c_s)))*(r(2)^2.d0);
T2_6 = (1.d0 + (r(2)/c_s))*(P1_bw - P_stat + P_v - P1_ext)/rho;
T2_8 = ( (-3.d0*poly_exp*r(2)*(P_stat - P_v + (2.d0*sigma/r1_eq))*((r1_eq/r(1))^(3.d0*poly_exp)) ) + (2.d0*sigma*r(2)/r(1)) - (4.d0*mu*(- ((r(2)^2.d0)/r(1)))) )/r(1);
T2_9 = 2.d0*pi*f_s*P_s*(cos(2.d0*pi*f_s*(t + (r(1)/c_s))))*(1.d0 + (r(2)/c_s) );
T2_7 = (r(1)/(rho*c_s))*(T2_8 - T2_9);
rdot(2) =  (T2_6 + T2_7 - T2_1 - T2_5)/(T2_4 + (4.d0*mu/(rho*c_s)));


rdot(3) = r(4);
P2_bw = ( (P_stat - P_v + (2.d0*sigma/r1_eq))*((r1_eq/r(3))^(3.d0*poly_exp)) ) - (2.d0*sigma/r(3)) - (4.d0*mu*r(4)/r(3));
P2_ext = P_s*sin(2.d0*pi*f_s*(t + (r(3)/c_s)));
T1_1 = ((2.d0*r(1)*(r(2)^2.d0)) + ((r(1)^2.d0)*rdot(2)))/d;
T1_4 = (1.d0 - (r(4)/c_s))*r(3);
T1_5 = 1.5d0*(1.d0 - (r(4)/(3.d0*c_s)))*(r(4)^2.d0);
T1_6 = (1.d0 + (r(4)/c_s))*(P2_bw - P_stat + P_v - P2_ext)/rho;
T1_8 = ( (-3.d0*poly_exp*r(4)*(P_stat - P_v + (2.d0*sigma/r1_eq))*((r1_eq/r(3))^(3.d0*poly_exp)) ) + (2.d0*sigma*r(4)/r(3)) - (4.d0*mu*(- ((r(4)^2.d0)/r(3)))) )/r(3);
T1_9 = 2.d0*pi*f_s*P_s*(cos(2.d0*pi*f_s*(t + (r(3)/c_s))))*(1.d0 + (r(4)/c_s) );
T1_7 = (r(3)/(rho*c_s))*(T1_8 - T1_9);
rdot(4) =  (T1_6 + T1_7 - T1_1 - T1_5)/(T1_4 + (4.d0*mu/(rho*c_s)));

rdot = rdot';

clc;
clear all;
close all;
time_range = [0 3000d-6];
initial_conditions = [10d-6 0.d0 5d-6 0.d0];
[t, r] = ode45('bubble', time_range, initial_conditions);
plot(t, r(:, 1), t, r(:, 3));

Ответы [ 3 ]

0 голосов
/ 25 октября 2018

По сути, у вас есть ситуация, когда 

rdot(2) = a2 + b2*rdot(4)
rdot(4) = a4 + b4*rdot(2)

, где a2,b2,a4,b4 содержит все другие термины в ваших выражениях.

Это линейная система, которую вам нужно решить, чтобы получитьправильные значения для возврата.Вы можете использовать линейный решатель Matlab или в этом простом двумерном случае сделать это вручную, 

rdot(2) = a2 + b2*(a4 + b4*rdot(2))  ==>  rdot(2) = (a2 + b2*a4) / (1 - b2*b4)
rdot(4) = a4 + b4*(a2 + b2*rdot(4))  ==>  rdot(2) = (a4 + b4*a2) / (1 - b2*b4)

Чтобы применить это, вам нужно разделить 

T1_1  as T1_1a + T1_1b*rdot(4),

Вы можете вычислить T1_1a и T1_1b из заданных констант и переменных состояния.Затем, где у вас есть в конце 

 rdot(2)=(other + coeff*T1_1)/denom

, вы должны разделить на 

 a2 =(other + coeff*T1_1a) / denom    and 
 b2 =         coeff*T1_1b  / denom

и сделать то же самое для второй части, чтобы получить a4,b4, а затем применить формулы решения выше,

0 голосов
/ 25 октября 2018

Я исправил проблему, введя две новые переменные с начальными условиями, но затем они обновляются по мере выполнения кода.Новый код здесь с двумя переменными: r2ddot и r1ddot; 

function rdot = f(t, r)

P_stat = 1.01325d5;
P_v = 2.3388d3;
mu = 1.002d-3;
sigma = 72.8d-3;
c_s = 1481d0;
poly_exp = 1.4d0;
rho = 998.2071d0;
f_s = 20d3;
P_s = 1.01325d5;
r1_eq = 4d-6;
r2_eq = 5d-6;
d = 50*(r1_eq + r2_eq);

r2ddot = 0;
r1ddot = 0;

rdot(1) = r(2);
P2_bw = ( (P_stat - P_v + (2.d0*sigma/r2_eq))*((r2_eq/r(3))^(3.d0*poly_exp)) ) - (2.d0*sigma/r(3)) - (4.d0*mu*r(4)/r(3));
P2_ext = P_s*sin(2.d0*pi*f_s*(t + (r(3)/c_s)));
T1_1 = ((2.d0*r(1)*(r(2)^2.d0)) + ((r(1)^2.d0)*r1ddot))/d;
T1_4 = (1.d0 - (r(4)/c_s))*r(3);
T1_5 = 1.5d0*(1.d0 - (r(4)/(3.d0*c_s)))*(r(4)^2.d0);
T1_6 = (1.d0 + (r(4)/c_s))*(P2_bw - P_stat + P_v - P2_ext)/rho;
T1_8 = ( (-3.d0*poly_exp*r(4)*(P_stat - P_v + (2.d0*sigma/r2_eq))*((r2_eq/r(3))^(3.d0*poly_exp)) ) + (2.d0*sigma*r(4)/r(3)) - (4.d0*mu*(- ((r(4)^2.d0)/r(3)))) )/r(3);
T1_9 = 2.d0*pi*f_s*P_s*(cos(2.d0*pi*f_s*(t + (r(3)/c_s))))*(1.d0 + (r(4)/c_s) );
T1_7 = (r(3)/(rho*c_s))*(T1_8 - T1_9);
rdot(2) =  (T1_6 + T1_7 - T1_1 - T1_5)/(T1_4 + (4.d0*mu/(rho*c_s))) ;
r2ddot = rdot(2);

rdot(3) = r(4);
P1_bw = ( (P_stat - P_v + (2.d0*sigma/r1_eq))*((r1_eq/r(1))^(3.d0*poly_exp)) ) - (2.d0*sigma/r(1)) - (4.d0*mu*r(2)/r(1));
P1_ext = P_s*sin(2.d0*pi*f_s*(t + (r(1)/c_s)));
T2_1 = ((2.d0*r(3)*(r(4)^2.d0)) + ((r(3)^2.d0)*r2ddot))/d;
T2_4 = (1.d0 - (r(2)/c_s))*r(1);
T2_5 = 1.5d0*(1.d0 - (r(2)/(3.d0*c_s)))*(r(2)^2.d0);
T2_6 = (1.d0 + (r(2)/c_s))*(P1_bw - P_stat + P_v - P1_ext)/rho;
T2_8 = ( (-3.d0*poly_exp*r(2)*(P_stat - P_v + (2.d0*sigma/r1_eq))*((r1_eq/r(1))^(3.d0*poly_exp)) ) + (2.d0*sigma*r(2)/r(1)) - (4.d0*mu*(- ((r(2)^2.d0)/r(1)))) )/r(1);
T2_9 = 2.d0*pi*f_s*P_s*(cos(2.d0*pi*f_s*(t + (r(1)/c_s))))*(1.d0 + (r(2)/c_s) );
T2_7 = (r(1)/(rho*c_s))*(T2_8 - T2_9);
rdot(4) =  (T2_6 + T2_7 - T2_1 - T2_5)/(T2_4 + (4.d0*mu/(rho*c_s)));
r1ddot = rdot(4);
rdot = rdot';


clc;
clear all;
close all;
time_range = [0 1d-3];
initial_conditions = [4d-6 0.d0 5d-6 0.d0];
[t, r] = ode45('bubble', time_range, initial_conditions);
plot(t, r(:, 1), t, r(:, 3));

Но я не могу получить желаемый результат.На самом деле я пытаюсь воспроизвести результаты из прилагаемого документа, см. Уравнение 7. введите описание ссылки здесь

enter image description here

ВторойНабор уравнений можно получить, поменяв местами индексы 1 и 2.

Важное примечание: в последнем члене уравнения 7 в статье Меттина есть опечатка, которую можно проверить с помощью проверки размеров различныхсрок.Правильный последний член можно увидеть из другой статьи https://journals.aps.org/pre/abstract/10.1103/PhysRevE.83.066313 См. Последний член в уравнении (1) ниже.Не обращайте внимания на другие дополнительные термины в уравнении.Важным моментом является то, что уравнение имеет второй порядок по Rj, а уравнение имеет член второго порядка по Ri в конце.И это то, что я пытался кодировать.enter image description here

Любая помощь будет высоко оценена.

0 голосов
/ 25 октября 2018

Для каждой степени ODE вам понадобится одно начальное условие.Это связано с количеством функций, которые вы рассчитываете.В вашем случае вы получили систему ODE второго порядка, которая после решения будет обеспечивать четыре функции: r(1), r(2), r(3), r(4)

Зачем нам даже нужны начальные условия?Представьте, что у вас есть простая производная: y'= y Мы знаем, что функция y = exp(x) * C решает эту проблему, но нам нужно настроить C, чтобы получить «Единственное решение».С другой стороны, нет смысла давать y' начальное условие, поскольку оно полностью определено после определения y.Не имеет значения, появляется ли «чужая» переменная в форме производной или в виде линейного фактора.Он не зависит от количества микросхем.

Я надеюсь, что смогу немного прояснить это.С моей точки зрения, ваша программа должна работать таким образом, но у меня не было возможности опробовать ее.

...