Решение диф уравнения с функцией, зависящей от решения

@zillant · Регистрация: 15.01.2014

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

Есть небольшая баллистическая программа. Осознал я недавно, что использовать постоянную плотность воздуха в атмосфере как минимум-глупо. Каких адекватных идей по реализации Pv(Y(:2)) у меня не появилось. Была идея с рекурсией, но что-то она не реализовалась. По сути. у меня получается какой-то замкнутый цикл, т.е. чтобы найти высоту, надо найти давление, а чтобы найти давление, надо найти высоту. Эмпирическая формула давления: P=2000-Y/2000+Y, где Y-высота.
Надеюсь, что кто-то из форумчанинов окажет помощь и подскажет)

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
function ballyst1
clear all
global Cx S Pv m a V0 g y
V0 = 780; % начальная скорость, м/с
a = 45; % угол возвышения орудия
g = 9.81; % ускорение свободного падения
m = 45; % масса снаряда, кг
D = 152; % калибр, мм
S = pi * ( D / 1000 ) ^ 2 % площадь поперечного сечения снаряда (миделево сечение)
Cx = 0.1; % аэродинамический(баллистический) коэффициент
Pv = 1.225; % плотность воздуха (на уровне моря)
NU = [0; 0]; % координаты в начальный момент времени
t = 100; % время, в течение которого ведется расчет
hold
%с сопротивлением
[T2,Y2]=ode45(@dV, [0 t], NU);
plot(Y2(:, 1), Y2(:, 2), '-g')
xlabel 'Range'
ylabel 'Height'
axis([0 50000 0 20000]); 
disp([T2,Y2]);
grid

Сама функция:

Matlab M
1
2
3
function dV = dV(t, dV)
global Cx S Pv  m a V0 g
dV = [(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];

@SSC · 05.12.2016, 07:10

Сообщение от zillant

Осознал я недавно, что использовать постоянную плотность воздуха в атмосфере как минимум-глупо.

Проблема не только в этом.
В функции dV при определении текущих ускорений по дальности и высоте Вы не используете текущих значений высоты и дальности, и кроме того при учете аэродинамического сопротивления расчитываете его не от текущей скорости.
А учет плотности воздуха - чего уж проще у Вас есть высота в текущий момент времени, используете Вашу формулу и получаете давление.
(А если у Вас нет высоты в текущий момент времени, то Вы списали программу и совершенно не понимаете что у Вас в ней написано)

@zillant · 05.12.2016, 12:45 **[ТС]**

Да, Вы правы, исользую не ту формулу, ошибся. Надо было сразу залезть в учебник Дмитриевского и не городить. Задачу эту я писал месяц назад где-то, и похоже что эту формулу выводил сам, правильнее будет делать как говорите Вы.
Изначально, я планировал, что записываю все параметры в массив, в формате (время, Х, У). В матлабе работать я начал относительно недавно, в основном использовал его как хороший калькулятор для лабораторных работ. В Дьяконове ничего похожего найти не смог, поэтому обратился сюда. Собственно, делал программу по примеру найденному в литературе.
Если использовать в программе обычный цикл, то как использовать текущее значение определенной величины мне понятно. А вот как его использовать при численном решении диф уравнения - не особо. Была идея передавать значение давления в начальных условиях, а потом его считать по предыдущей высоте, но опять же, реализация хромает.

@SSC · 05.12.2016, 13:35

Сообщение от zillant

function dV = dV(t, dV)

Использование одного имени переменной сразу для возвращаемого значения, имени функции и входного параметра - это несколько перебор.
Это как двух братьев близнецов назвать Вася, а потом еще и их сестру Васелисой, и потом по фразе Вася П... не сделал уроки, Вы сразу точно знаете кто виновник торжества.
Если обозначить для примера так

Matlab M
1
function dV = dVf(t, z)

тогда z(0) и z(1) будут дальностью и высотой (или наоборот) в расчетный момент времени.
Кроме этого желательно перейти на диффуравнения второго порядка для дальности и высоты - это упростит расчеты.
Посмотри MATLAB-овский пример прыгающего мяча http://www.mathworks.com/help/... ation.html

Новые блоги и статьи Все статьи Все блоги /
делаю науч статью по влиянию грибов на сукцессию anaschu 13.03.2026 прикрепляю статью	SDL3 для Desktop (MinGW): Создаём пустое окно с нуля для 2D-графики на SDL3, Си и C++ 8Observer8 10.03.2026 Содержание блога Финальные проекты на Си и на C++: hello-sdl3-c. zip hello-sdl3-cpp. zip Результат:	Установка CMake и MinGW 13.1 для сборки С и C++ приложений из консоли и из Qt Creator в EXE 8Observer8 10.03.2026 Содержание блога MinGW - это коллекция инструментов для сборки приложений в EXE. CMake - это система сборки приложений. Здесь описаны базовые шаги для старта программирования с помощью CMake и. . .	Как дизайн сайта влияет на конверсию: 7 решений, которые реально повышают заявки Neotwalker 08.03.2026 Многие до сих пор воспринимают дизайн сайта как “красивую оболочку”. На практике всё иначе: дизайн напрямую влияет на то, оставит человек заявку или уйдёт через несколько секунд. Даже если у вас. . .
Модульная разработка через nuget packages DevAlt 07.03.2026 Сложившийся в . Net-среде способ разработки чаще всего предполагает монорепозиторий в котором находятся все исходники. При создании нового решения, мы просто добавляем нужные проекты и имеем. . .	Модульный подход на примере F# DevAlt 06.03.2026 В блоге дяди Боба наткнулся на такое определение: В этой книге («Подход, основанный на вариантах использования») Ивар утверждает, что архитектура программного обеспечения — это структуры,. . .	Управление камерой с помощью скрипта OrbitControls.js на Three.js: Вращение, зум и панорамирование 8Observer8 05.03.2026 Содержание блога Финальная демка в браузере работает на Desktop и мобильных браузерах. Итоговый код: orbit-controls-threejs-js. zip. Сканируйте QR-код на мобильном. Вращайте камеру одним пальцем,. . .	SDL3 для Web (WebAssembly): Синхронизация спрайтов SDL3 и тел Box2D 8Observer8 04.03.2026 Содержание блога Финальная демка в браузере. Итоговый код: finish-sync-physics-sprites-sdl3-c. zip На первой гифке отладочные линии отключены, а на второй включены:. . .

@zillant 0 / 0 / 0 Регистрация: 15.01.2014 Сообщений: 49
	05.12.2016, 12:45 [ТС]
	Да, Вы правы, исользую не ту формулу, ошибся. Надо было сразу залезть в учебник Дмитриевского и не городить. Задачу эту я писал месяц назад где-то, и похоже что эту формулу выводил сам, правильнее будет делать как говорите Вы. Изначально, я планировал, что записываю все параметры в массив, в формате (время, Х, У). В матлабе работать я начал относительно недавно, в основном использовал его как хороший калькулятор для лабораторных работ. В Дьяконове ничего похожего найти не смог, поэтому обратился сюда. Собственно, делал программу по примеру найденному в литературе. Если использовать в программе обычный цикл, то как использовать текущее значение определенной величины мне понятно. А вот как его использовать при численном решении диф уравнения - не особо. Была идея передавать значение давления в начальных условиях, а потом его считать по предыдущей высоте, но опять же, реализация хромает. 0