MATLAB решение дифура

[zillant]

Есть такая система, если учесть что плотность P постоянная величина, то все получается чудесно, но если считать его по барометрической формуле, то функция отказывается работать. Как мне использовать метод решения диф уравнения, чтобы там была переменная величина. Основная трудность заключается в том, что плотность зависит от высоты, а я ее как раз и ищу методом решения диф уравнения.

Используется синтаксис Matlab M

[T2,Y2]=ode45(@Complex, [0 t], NU);
function Complex = Complex(t, x)
global Cx S Pv m a V0 g
Pv = 0.0029*101.325*exp(-0.00012*y)/8.31*293;
Complex = [ ( 2*m*V0*cosd(a) ) / ( Cx*S*Pv*V0*cosd(a).*t + 2*m ) ; ( 2*m*V0*sind(a) ) / ( Cx*S*Pv*V0*sind(a).*t + 2*m ) - g.*t ];

В массиве Complex две составляющих, Vx и Vy

[Pphantom]

У Вас Pv зависит от y, который больше нигде и никак не упоминается. Откуда он должен взяться?

[zillant]

[T2,Y2]=ode45(@Complex, [0 t], NU);
в свою очередь Y2=[X,Y]

Я как раз его и нахожу, решая диф уравнение. Можно ли как решить такое диф уравнение, одна из величин которого зависит от решения.

[Pphantom]

То, что Вы назвали x в функции, вообще говоря, должно представлять собой вектор-функцию, одним из компонентов которой является то, что названо y.

[zillant]

то что я передаю то в функцию, да. А вот как мне сделать так чтобы туда передавать туда как в цикле значения, и уже исходя от этого диффур решался?

[Pphantom]

zillant в сообщении #1181892 писал(а):

то что я передаю то в функцию, да. А вот как мне сделать так чтобы туда передавать туда как в цикле значения, и уже исходя от этого диффур решался?

Зачем? Вы используете уже готовый интегратор, не надо пытаться делать за него его работу.

[zillant]

Тогда как мне быть? Если делать два диффура, в одном Pv постоянная, получить там значения для Y, а потом передавать их в диффур с Pv(Y), то мне кажется с точностью будет не очень хорошо...

[Pphantom]

Может быть, просто почитать справку по ode45? Ну или для начала замените y в заголовке функции на x(2), если, конечно, начальные условия и т.п. Вы задали правильно.

[zillant]

Не поверите, много чего изучил, проштудировал Дьяконова, множество информации прогуглил, и прише к отчаянию, задал вопрос на форуме

function Complex = Complex(t, x)
global Cx S  m a V0 g
%Pv = 0.0029*101.325*exp(-0.00012*x(2))/8.31*293;
Complex = [ ( 2*m*V0*cosd(a) ) / ( Cx*S*0.0029*101.325*exp(-0.00012*x(2))/8.31*293*V0*cosd(a).*t + 2*m ) ; ( 2*m*V0*sind(a) ) / ( Cx*S*0.0029*101.325*exp(-0.00012*x(2))/8.31*293*V0*sind(a).*t + 2*m ) - g.*t ];

Результаты получаются вообще не очень, работает не так как надо.

[zillant]

Нужно построить траекторию движения снаряда, с учетом сопротивления воздуха. Плотность воздуха считать зависимой от высоты.
Я вывел следующие формулы:
$Pv=0.0029\cdot101.325\exp(-0.00012Y)/8.31\cdot293$
Vx=$\frac{2mV0cos(a)}{CSPvV0cos(a)t + 2m }$
Vy=$\frac{2mV0sin(a)}{CSPvV0sin(a)t + 2m }-gt$
Теперь нужно решить систему дифуравнений для получения координат.

Используется синтаксис Matlab M

V0 = 680; % начальная скорость
a = 55; % угол бросания (возвышения орудия)
g = 9.8; % ускорение свободного падения
m = 10; % масса снаряда
cal = 76; % калибр (в мм)
S = pi * ( cal / 1000 ) ^ 2 % площадь поперечного сечения снаряда (мидель)
C = 0.1; % аэродинамический(баллистический) коэффициент

[GAA]

У Вас изначально в координатной записи система двух уравнений второго порядка. Уравнения содержат вторые производные координат по времени ($d^2x/dt^2$, $d^2y/dt^2$), одну координату ($y$) и проекции скорости ($dx/dt$, $dy/dt$) движения снаряда. Для использования ode45, задача должна быть сведена к системе четырех уравнений первого порядка. Не ясно как Вы смогли её свести к двум уравнениям первого порядка. Лучше все же записать исходную задачу. И указать каким образом Вы её решали.

[zillant]

Как таковой исходной задачи нет. Хочу построить траекторию и все. Да, я вручную рассмотрел два уравнения второго порядка. Получаются из 2ЗН.
$m\frac{dV_x}{dt}=\frac{-P_vV_x^2SC}{2}$
$m\frac{dV_y}{dt}=\frac{-P_vV_y^2SC}{2}-mg$
После, от них вручную интегрирую и получаю следующие уравнения(интегрирую со знанием того, что здесь вторая и первая производная). Изначально делал задачу с постоянной плотностью воздуха на уровне моря, все работало отлично. Сейчас я хочу усовершенствовать модель, добавив плотность воздуха зависящей от высоты. Формулу плотности я вывел из барометрической формулы.
Возможно вы не так поняли, при записи формул я допустил небольшую оплошность:
$P_v=0.0029\cdot101.325\exp(-0.00012Y)/8.31\cdot293$
$V_x=\frac{2mV_0cos(a)}{CSP_vV_0cos(a)t + 2m }$
$V_y=\frac{2mV_0sin(a)}{CSP_vV_0sin(a)t + 2m }-gt$

[GAA]

Ваша система должна содержать ещё два уравнения:
$dx/dt = V_x$
$dy/dt = V_y$
Мы не можем проинтегрировать отдельно два выписанных в предыдущем Вашем сообщении уравнения:
$m\frac{dV_x}{dt}=\frac{-P_vV_x^2SC}{2}$
$m\frac{dV_y}{dt}=\frac{-P_vV_y^2SC}{2}-mg$

-- Чт 05.01.2017 11:50:05 --

Т.е. уравнения для компонент скорости не содержат независимую переменную. Можно выполнить стандартную замену и, используя выражение для $P_v$, записать квадратуры (не пробовал брать эти квадратуры, возможны затруднения). Но ожидаемый результат не будет похож на Ваши выражения для скоростей: результат интегрирования будет связывать скорости и координаты, но не будет содержать время. (Зависимость от времени при таком подходе появится при последующем интегрировании уравнений $dx/dt = V_x$, $dy/dt = V_y$.)

Но я сразу бы считал в Matlab систему четырех уравнений. Если нужно только как-то построить траектории.

[zillant]

Вы предлагаете что-то вроде этого?

Используется синтаксис Matlab M

function dVdt = trajecroty (t,coords)
dVdt(1,1)=[x(2);y(2)]
dVdt(2,1) = [((-0.0029*101.325*exp(-0.00012*coords(2)/8.31*293)*(x(2))^2*S*C)/(2*m);(-(0.0029*101.325*exp(-0.00012*coords(2))/8.31*293)*(y(2))^2*S*C)/(2*m)-g];
 

[GAA]

По приведенному Вами тексту я ничего не могу понять. Но сведение двух обыкновенных дифференциальных уравнений второго порядка к системе четырех уравнений первого порядка - стандартный материал. (См. любой учебник, например, Эльсгольц "Обыкновенные дифференциальные уравнения и вариационное исчисление".) Если вопрос по дифференциальным уравнениям, то набирайте формулы в $\TeX$.

Если Вы спрашиваете о реализации в Matlab, то фрагмент слишком короткий, чтобы что-то содержательное говорить. (Но такое я ни разу не видел.) Приводите более объёмные исходные тексты. Посмотреть можно, например, в примере 14 (движение тела, брошенного с начальной скоростью $v_0$ под углом $\alpha$ к горизонту в предположении, что сопротивление воздуха пропорционально квадрату скорости) Доля П.Г. "Использование MATLAB. Решение дифференциальных уравнений." (В Сети был доступен pdf).