Scilab в свободном падении / forpes.ru

Главная
Scilab в свободном падении

Scilab в свободном падении +42

22.04.2017 11:47

shbma 62 7200 Источник

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

Как известно, свободное падение — это равнопеременное движение под действием силы тяжести, когда другие силы, действующие на тело отсутствуют или пренебрежимо малы. В быту часто свободным падением также называют движение в атмосфере Земли, когда на ускорение тела не влияют никакие силы, кроме силы сопротивления воздуха и силы тяжести. Именно таким расширенным пониманием мы и будем пользоваться в дальнейшем.

Задача

Составить диф.уравнения (ДУ), описывающие существенные свойства полета Баумгартнера, решить их средствами Scilab и получить в явном виде зависимости дистанции и скорости от времени свободного падения.

Исходные данные

Т.к. основная задача статьи — все-таки показать, как Scilab умеет решать ОДУ, а не получить идеально точную модель, увлекаться анализом экспериментальных данных (телеметрии), погрешностей и подгонкой мы особенно не будем. Однако для проверки адекватности модели несколько значений все же необходимо взять. Итак:

Высота прыжка — 39000м.
Дистанция свободного падения — 36402.6м.
Время свободного падения — 4 мин 20 сек.
За первые 20сек Баумгартнер набрал скорость 194 м/сек, через 50сек падения (на высоте 8447м) он развил рекордную скорость в 377м/сек, а к моменту раскрытия парашюта (т.е. через 260 сек полета) скорость падения снизилась до 77 м/сек.

Источники: Википедия, Пресса, Видео.

Модель-1 (без учета сопротивления воздуха)

Начнем с простейшей модели, в которой есть только одна сила — сила притяжения Земли. Полагаем, что падение идет строго вертикально, начало координат размещено в точке начала прыжка, ось-y направлена вниз.

Ускорение экстремала, соответственно, равно ускорению свободного падения.

$\frac{dv}{dt} = g$

При этом по определению скорости:

$\frac{dy}{dt} = v$

Получили систему из двух обыкновенных ДУ. Решаем.

Воспользуемся функцией ode(), входящей в стандартный дистрибутив Scilab-а.
Согласно документации, ode решает явные обыкновенные дифференциальные уравнения, определенные как:

$\frac{dy}{dt} =f(t,y)\\ y(t_0) = y_0$

Прямо как в нашем случае. Слева должны быть производные 1го порядка.

Если же ДУ содержит производную 2го, 3го и пр. порядков, то нужно ввести замену(ы) и преобразовать одно уравнение с систему нескольких. Как мы на самом деле и сделали. Ведь ускорение — это 2я производная координаты по времени, от которой мы перешли к 1й производной, введя замену — переменную «скорость» равную dy/dt. Т.е.
от $\frac{d^2y}{dt^2} = g$ перешли к $$inline$\begin{equation*} \begin{cases} \frac{dv}{dt} = g \\ \frac{dy}{dt} = v \end{cases} \end{equation*}$inline$$

Наконец-то переходим к коду (ссылка на github). Его можно писать в предлагаемом средой редакторе SciNotes, равно как и в любом другом текстовом редакторе. Можно также запустить оболочку без графики и вводить команды в консоли. На мой взгляд SciNotes удобнее подсветкой кода и интеграцией в среду разработки.

Первым делом нам нужно описать ДУ в виде функции на языке Scilab.

//Функция, описывающая систему ДУ динамики падающего в вакууме тела
//vector = [v,y], t - время
function dVector = verticalFallingVacuum(t, vector)
    dVector = zeros(2,1);   //заготовка под результат - вектор из двух  нулей
    
    dVector(1) = g;         //dv/dt = g
    dVector(2) = vector(1); //dy/dt = v
endfunction

Затем задаем начальные условия.

//Граничные и начальные условия
y0 = [0;0];     //v и x в начале движения, м/сек и сек.
t0 = 0;         //начальный момент времени
t = 0:0.1:260;  //временной интервал движения, сек. (4мин 20сек)

После чего вызываем функцию-решатель ode(), передав ей на вход все созданное выше.

result = ode(y0,t0,t,verticalFallingVacuum); //решение ДУ

В матрице result в итоге окажуться значения искомых функций-решений данной системы ДУ, соответствующие точкам оси времени из указанного в t диапазона.

Посмотрим на то, что получилось в графической форме, а заодно выведем в консоль пару контрольных точек.

//строим График изменения пройденной дистанции от времени
plot(t,result(2,:)); //вторая строка матрицы, по всему диапазону
xtitle("Дистанция свободного падения как функция от Времени"); //заголовок графика
//выведем в консоль пару контрольных значений
mprintf("Cкорость через 50сек = %f км/ч \n",result(1,500)*3.6);
mprintf("Пройденный путь за 4мин.20сек = %f км \n",result(2,2600)/1000);

Все вместе:

g = 9.81; // ускорение свободного падения м/с

mode(0);    //разрешаем вывод в консоль

//Функция, описывающая систему ДУ динамики падающего в вакууме тела
//vector = [v,y]
function dVector=verticalFallingVacuum(t, vector)
    dVector = zeros(2,1);   //заготовка под результат
    
    dVector(1) = g;         //dv/dt = g
    dVector(2) = vector(1); //dy/dt = v
endfunction

//Граничные и начальные условия
y0 = [0;0];     //v и x в начале движения, м/сек и сек.
t0 = 0;         //начальный момент времени
t = 0:0.1:260;  //временной интервал движения, сек. (4мин 20сек)
result = ode(y0,t0,t,verticalFallingVacuum); //решение ДУ

//строим График
plot(t,result(2,:));
xtitle("Дистанция свободного падения как функция от Времени")
//выведем в консоль пару контрольных значений
mprintf("Cкорость через 50сек = %f км/ч \n",result(1,500)*3.6);
mprintf("Пройденный путь за 4мин.20сек = %f км \n",result(2,2600)/1000);

Получилась ожидаемая квадратичная зависимость координаты от времени.

Сверимся с практикой:
Cкорость через 50сек = 1762 км/ч (а должно быть 1357 км/ч).
Пройденный путь за 4мин.20сек = 331.3 км (а должно быть 36.4 км.).

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

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

Функция, строящая улучшенный график

Функция-графопостроитель:

//Строит графики зависимостей скорости и координаты от времени свободного падения.
//t - массив 1xN точек на временной оси
//data - расчетные данные - массив 2xN точек v,y
//realJumpVT - массив 2хN - экспериментальные точки скорость-время.
//realJumpYT - массив 2xN экспериментальные точки координата-время.
function [yAxis,vAxis] = drawGraphics (t, data, realJumpVT, realJumpYT)
    scf();//создаем графическое окно

    xsetech([0 0 0.9 0.47]);//создаем под-окно - выделяем верхнюю половину окна под график
    plot2d(t', data(2,:)'); 
    
    yAxis=gca(); // возьмем дескриптор текущего гарфика
    yAxis.title.text = "Дистанция свободного падения, м";
    yAxis.title.font_size=2; //увеличим размер шрифта заголовка
    yAxis.x_label.text = "Время с момента прыжка, с";
    yAxis.y_label.text = "Дистанция от точки прыжка, м";
    yAxis.grid = [0,0];    //включаем сетку и длаем ее черной
    yAxis.x_location = "origin"; // пустим x-ось через ноль
    yAxis.y_location = "origin"; // пустим y-ось через ноль
    yAxis.children(1).children(1).foreground = 11; //синий цвет линии
    
    plot2d(realJumpYT(2,:), realJumpYT(1,:)',[-2]); //наносим экспериментальные точки
    grExpY=gca(); // возьмем дескриптор текущего гарфика
    grExpY.children(1).children(1).mark_foreground = 14; //маркер - зеленый
    grExpY.children(1).children(1).thickness = 2; 
    
    leg1 = legend(["теория", "реальная жизнь"],-1);
    
    xsetech([0 0.53 0.9 0.43]);//создаем под-окно - выделяем нижнюю половину окна под график
    plot2d(t', data(1,:)'); 
    vAxis=gca(); // возьмем дескриптор текущего гарфика
    vAxis.title.text = "Скорость свободного падения, м/с";
    vAxis.title.font_size=2; //увеличим размер шрифта заголовка
    vAxis.x_label.text = "Время с момента прыжка, с";
    vAxis.y_label.text = "Скорость падения, м/с";
    vAxis.grid = [0,0];    //включаем сетку и длаем ее черной
    vAxis.x_location = "origin"; // пустим x-ось через ноль
    vAxis.y_location = "origin"; // пустим y-ось через ноль
    vAxis.children(1).children(1).foreground = 11; //линяя - синяя
    
    plot2d(realJumpVT(2,:), realJumpVT(1,:)',[-2]); //наносим экспериментальные точки
    grExpV=gca(); // возьмем дескриптор текущего гарфика
    grExpV.children(1).children(1).mark_foreground = 14; //маркер - зеленый
    grExpV.children(1).children(1).thickness = 2; 
    
    plot2d(t', 300*ones(size(t,'c'),1),[3]); //рисуем звуковой барьер
    grBarrier=gca(); // возьмем дескриптор текущего гарфика
    grBarrier.children(1).children(1).foreground = 5; //линия - красная
    xstring(200,320,"Звуковой барьер");
    t=get("hdl")   //получаем дискриптор только что созданного объекта
    //t.font_color=21; // бордовый
    t.font_size=2;
    
    leg2 = legend(["теория", "реальная жизнь"],-1);

endfunction

Исправим ситуацию, введя в модель силу сопротивления воздуха.

Модель-2 (с учетом сопротивления воздуха)

Атмосфера будет действовать на летящего с силой

$F_{air} = c \cdot S \cdot density \cdot \frac{v^2}{2},$

где
с — коэффициент сопротивления
S — наибольшая площадь сечения, перпендикулярного направлению полета
density — плотность воздуха
v — скорость движения
Формула взята отсюда, где также можно почитать аналитическое решение

Нас интересует ускорение (хотя лучше было бы сказать замедление), которое F_air будет сообщать летящему вниз экстремалу, чтобы добавить поправку в первое уравнение системы.

По 2му закону Ньютона $F_{air} = m \cdot a$

$a = \frac{F_{air}}{m} = \frac{c \cdot S \cdot deinsity \cdot v^2}{2 \cdot m}$

С учетом данной поправки наша система примет вид

$$$display$$\begin{equation*} \begin{cases} \frac{dv}{dt} = g - \frac{1}{m} \cdot \frac{c \cdot S \cdot deinsity \cdot v^2}{2}\\ \frac{dy}{dt} = v \end{cases} \end{equation*}$$display$$$

Решаем аналогичным образом. Описываем систему ДУ как Scilab-функцию.

//Функция, описывающая систему ДУ динамики падающего тела с учетом сопротивления воздуха 
//vector = [v,y], t-время
function dVector = verticalFallingAir(t, vector)   
    dVector = zeros(2,1);   //вектор-болванка под будущий результат
    
    height = startHeight-vector(2); //переводим дистанцию св.падения в высоту от земли
    aAir = (0.5*c*airDensity(height)*S*vector(1)^2)/m; //"ускорение" за счет торможения о воздух
    dVector(1) = g - aAtm;      //dv/dt = g - aAtm
    dVector(2) = vector(1);     //dy/dt = v    
endfunction

Поскольку плотность атмосферы изменяется в зависимости от высоты, нам потребуется написать соответствующую функцию. Таблицу можно взять отсюда. А вообще есть ГОСТ 4401-81 Атмосфера стандартная. Параметры (с Изменением N 1).

//Вычисляет плотность атмосферы на заданной высоте
function density = airDensity(height)
    hVSdensity = [
        0         1.225;
        50       1.219;
        100     1.213;
        200     1.202;
        300     1.190;
        500     1.167;
        1000     1.112;
        2000    1.007;
        3000    0.909;
        4000    0.819;
        5000    0.736;
        8000    0.526;
        10000   0.414;
        12000   0.312;
        15000   0.195;
        20000   0.089;
        50000   1.027*10^(-3);
        100000  5.550*10^(-7);
        120000  2.440*10^(-8);
    ];
    //проверим высоту
    if (height < 0) then
        height = 0;
    end
    //линейно интерполируем и берем значение на заданной высоте
    density = interpln(hVSdensity',[height]); 
endfunction

Задаем начальные условия и решаем.

//Граничные и начальные условия
step = 0.1;
y0 = [0;0];       //v и x в начале движения, м/сек и сек.
t0 = 0;           //начальный момент времени
t = 0:step:260;   //временной интервал движения, сек. (4мин 20сек)

result = ode(y0,t0,t,verticalFallingAir); //решение ДУ

Построим график и посмотрим, что получилось.

Видим, что между 40-й и 70-й секундами пробьем звуковой барьер, в районе 50й секунды скорость достигнет максимума, после чего начнет снижаться. Набранной скорости не хватит, чтобы сгореть в атмосфере (на Хабре обсуждалось), кроме того к концу полета скорость снизится до величины, позволяющей безопасно открыть парашют. Что соответствует реальности.

Скорость — сравнение теории с экспериментом в контрольных точках

Время св.падения, с	Расчетная скорость, м/с	Фактическая скорость, м/с
20	183	194
50	312	377
260	73	77

Дистанция св.падения — сравнение теории с экспериментом в контрольных точках

Время св.падения, с	Расчетная дистанция, м	Фактическая дистанция, м
50	9853	8447
260	41049	36403

По скорости ошибка в среднем 27 м/с (7% от фактического рекорда), по дистанции свободного падения 3026 м (8% от фактически пройденной в св.падении дистанции). Для наших учебно-демонстрационных целей вполне приемлемо.

Заключение

Таким образом, мы на примере простой, но наглядной, мат. модели рассмотрели как при помощи Scilab решать ДУ, интерполировать, строить графики. Надеюсь эта статься будет полезна начинающим изучать Scilab, а кого-то быть может и подвигнет на дальнейшее применение Scilab в своей практической деятельности.

Поделиться с друзьями

-->

Комментарии (62)

lzb_j77
22.04.2017 15:29
#10186198
Пытался я использовать Scilab — падает постоянно. Пришлось удалить.
1. Kazancev
  23.04.2017 14:51
  #10187254
  Наблюдал падения на версиях ниже 6-й беты при выделении памяти больше, чем задавалось в переменной стека (stacksize), и скорее всего глюки были связаны именно с ограниченным стеком, ну и java-машина что-то накручивала. В 6-й, пока, не замечал, так как память не переполнял — она теперь неограничена по-сути.
  1. CrazyFizik
    24.04.2017 10:51
    #10188174
    Кстати да. А что там сейчас действительно с памятью?
    Т.к. я тащусь от больших матриц, а вот Scilab (по крайней мере раньше) — нет
    
    Kazancev
    24.04.2017 18:44
    #10189038
    Смотря что Вы понимаете под большими матрицами. В версиях ниже 6-й четыре полные матрицы 8000х8000 забивали всю доступную память — выделялось ~ 400 МБ максимум если не ошибаюсь, — а уже 6-я бета проблем с памятью не имела, как минимум доступна вся RAM (файл подкачки не проверял). Тем не менее мне выгоднее использовать разрежённые матрицы, sparse, — я привёл скриншот в конце всех комментариев к этой статье, где формируются sparse-массивы.
    
    CrazyFizik
    24.04.2017 22:52
    #10189380
    Хз… ну там матрица ранга эдак 100к :-)
    Ну или хотя бы 15к
    
    Еще такой вопрос. А как там с поддержкой распараллеливания и CUDA?
    
    Kazancev
    25.04.2017 21:15
    #10191010
    Есть модуль GPU, не использовал. Возможно будет проще написать код на Fortran/Си, создать библиотеку и подцепить в Scilab. К слову, однопоточный битарник от Fortran'а с максимальной оптимизацией выполняется только в 1.5-2 раза быстрее, чем код в Scilab. Выполняются ли стандартные функции Scilab в многопотоке — не в курсе.
    
    Версия
    2.0-0
    Авторы
    Cedric Delamarre Vincent LEJEUNE, Allan CORNET, Cedric DELAMARRE
    Описание
    This module provides GPU computing capabilities to Scilab. It uses an implementation of BLAS (cuBLAS) and FFT (cuFFT) through gpuAdd, gpuMult, gpuFFT and other functions. This module uses essentially Cuda but some functions, as gpuBuild, have been created for build and use kernels developed with OpenCL or Cuda. Binary versions have been built with Cuda 4.0. If you need the OpenCL version of this module, please send an email to contact@scilab-enterprises.com.
    Смотри также
    • http://atoms.scilab.org/toolboxes/sciGPGPU/2.0 • http://forge.scilab.org/index.php/p/sciCuda/
    Дата выпуска
    2013-09-11
    Размер загр
1. vidyacat
  23.04.2017 22:30
  #10187738
  -2
  как мне сказали однажды: «это опенсорс, если мешают баги фикси их сам». поэтому лучше не используйте опенсорс.
  1. Kazancev
    24.04.2017 18:49
    #10189042
    Как только потребности пользователя перешагнут возможности OS и на горизонте замаячат убытки (недополученная прибыль) от его использования — пользователь купит коммерческий софт. «Экономика должна быть экономной». Сегодня мои потребности Scilab удовлетворяет, а что будет завтра — будет завтра.

SvyatoslavMC
22.04.2017 16:34
#10186244
+2
> В ряде ВУЗов (к примеру, УрФУ, ИТМО) ее используют для обучения студентов.

Почти во всех ВУЗах её используют, т.к. сейчас строже стали относится к обучению на софте, который студент технически может только спиратить откуда-нибудь. Сам застал переписывание методичек с Matlab на Scilab, и с MS Word на Libreoffice.

Scilab отличный инструмент и не падает постоянно, как у комментатора выше. Хотя некоторые падения могу подтвердить: в коде можно допустить ошибки, из-за который отвалится скрипт вместе со всей средой выполнения.
1. apro
  22.04.2017 18:37
  #10186310
  +3
  Сам застал переписывание методичек с Matlab на Scilab
  А почему не octave, у него вроде язык совместим с matlab,
  т.е. методички и переписывать бы не пришлось?
  1. SvyatoslavMC
    22.04.2017 19:20
    #10186334
    Я тогда не знал про Octave, а преподаватели может и до сих пор не знают.
1. lzb_j77
  23.04.2017 10:39
  #10186946
  Я скрипиты не писал сам, я брал примеры — падает, и всё тут. Нам же надо не глюки выискивать, а работать.
  Пришлось добывать Матлаб — там ничего не падает.
  1. SvyatoslavMC
    23.04.2017 11:01
    #10186974
    Там есть проблема совместимости. Написанные ранее скрипты могут просто падать в свежих версиях программы. Видимо, именно с этим вы столкнулись, когда брали чьи-то примеры.
    
    kasyachitche
    23.04.2017 11:28
    #10187010
    Там есть проблема совместимости. Написанные ранее скрипты могут просто падать в свежих версиях программы.
    
    С этой фразы следовало начинать пост.
1. altai2013
  23.04.2017 12:25
  #10187090
  +1
  Поделитесь, как человек, знакомый с этой вузовской кухней, почему ВУЗ не может просто купить хотя бы 20 лицензий в аудиторию? Ведь у разработчиков софта есть специальные льготные программы для учебных заведений, некоторые даже бесплатно дают им лицензии. Я понимаю техническую и экономическую сложность приобретения каких-нибудь банковских систем, но Matlab стоит всего $550 на пользователя и при оптовой покупке наверняка есть еще скидки. У ректоров зарплата по 600 тысяч в месяц и по миллиону, неужели нельзя раз в 10 лет приобрести пачку лицензий на базовый учебный софт?
  1. SvyatoslavMC
    23.04.2017 12:39
    #10187114
    +2
    Не знаю, какие у кого зарплаты, но на подобные предложения обычно ответ в силе «Денег нет, но вы держитесь».
    
    Ещё я пытался воспользоваться программой Microsoft DreamSpark, оплаченной ВУЗом. Тогда я потратил 3 месяца на поиск ответственных людей, но в итоге мне выдали несколько ключей, чтобы я отстал, хотя по документации, которую я даже распечатал и носил по кабинетам, должны были просто подтвердить мою учётную запись в админском аккаунте.

Ilyato
22.04.2017 18:35
#10186306
Любопытно. «Фактическая скорость» во всех 3х точках выше расчетной, а «фактическая дистанция» во второй и третей точках меньше расчетной.

shbma
22.04.2017 19:04
#10186318
У нас на Scilab переходили с MathCAD, так что о совместимости с MatLab никто особо не думал.

LongNowhere
22.04.2017 19:41
#10186360
+1
У нас тоже переходили с Matlab на Scilab.
Отсутствие анонимных функций утомляет после Matlab, а так весьма достойная свободная альтернатива.

oburdelev
22.04.2017 19:44
#10186368
-1
> g = 9.81;
Стал читать статью именно в ожидании увидеть интересное решение для интегрирования уравнения свободного падения. Не увидел. Хотя для этой задачи зависимость ускорения свободного падения от высоты даст значимый вклад.
1. nerudo
  22.04.2017 19:59
  #10186398
  При высоте 30 км разница g вроде не более 1% в верхней точке.

maisvendoo
22.04.2017 20:33
#10186450
+4
За статью плюс, но позволю себе замечание чисто методического характера.

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

$m \, \vec a = \sum \vec F_k$

произведение массы на ускорение есть векторная сумма действующих на точку сил. Так что Ваши выкладки, связанные с вычислением ускорения силы сопротивления воздуха избыточны. Достаточно написать

$m \, \vec a = m \, \vec g + \vec R$

и проецировать уравнение на вертикальную ось, сведя все выкладки к делению обоих частей получившегося уравнения на массу.

К тому же слово «замедление» это инженерная вульгарность. Ускорение величина векторная и знака у него нет. Знак есть у проекции, но проекция это не есть ускорение в механическом смысле
1. San_tit
  23.04.2017 13:02
  #10187134
  Производная всё-таки имеет знак, а замедление относится все-таки скорее к знаку второй производной, а не к вектору как таковому
  1. maisvendoo
    23.04.2017 13:47
    #10187176
    Производная всё-таки имеет знак, а замедление относится все-таки скорее к знаку второй производной, а не к вектору как таковому
    
    Производная от вектора по времени — тоже вектор, и никакого знака у ней нет. Не путайте понятие векторной величины и её проекции, производной от вектора, и производной от его проекции

progman_rus
22.04.2017 21:19
#10186502
а можно так рассчитать конечную скорость тела падающего на звезду массой М?
1. kasyachitche
  23.04.2017 11:34
  #10187020
  Для этого понадобится не один десяток параметров, да и вряд ли кто-то изучал требуемые свойства гелиосферы с достаточной точностью.
  
  А если что-то и посчитать, то будут ли основания верить полученному решению?
  1. progman_rus
    23.04.2017 16:41
    #10187342
    Если честно меня всегда интеерсовал один вопрос — какой максимальной скорости может достичь тело при падении на поверхность планеты/звезды/ЧД.
    Без учетов давления света и прочих факторов.
    Своих сил не хватает решить диффур — там g переменная от R (((
    
    kasyachitche
    23.04.2017 18:46
    #10187508
    Если проинтегрировать данное выражение (из Википедии)
    
    по h от -inf до 0, то получится:
    , где G — гравитационная постоянная, M — масса планеты, r — радиус планеты.
    Вот вам и максимальная скорость, при падении из максимально удаленной от планеты точки, до ее поверхности без учета сопротивления атмосферы и релятивистских эффектов.
    
    progman_rus
    23.04.2017 19:18
    #10187528
    о! пасибо!
    
    maisvendoo
    23.04.2017 19:31
    #10187538
    Всё уже проинтегрировали за нас ). Пусть точка падает на планету/звезду из «бесконечности» без начальной скорости. Потенциальная энергия точки в гравитационном поле равна
    
    $\Pi = - \frac{G\,mM}{r}$
    r — расстояние до притягивающего центра. За ноль потенциальной энергии принимаем положение точки на «бесконечности». В отсутствии сил сопротивления сохраняется полная механическая энергия точки, тогда
    
    $0 = \frac{m \, v^2}{2} - \frac{G\,mM}{R}$
    где v — скорость точки при падении на поверхность; R — радиус небесного тела. Получаем искомую скорость
    
    $v = \sqrt{\frac{2\,G\,M}{R}}$
    и это внезапно — вторая космическая скорость на поверхности небесного тела
    
    progman_rus
    23.04.2017 19:48
    #10187560
    интересно, а падая на очень массивную зведу можно превысить с
    ОТО запрещает… но по формулам получается что конечная скорость зависит только отмассы звезды и ее радиуса.
    интересно. наверное какой нибудь казус есть в ОТО с поправка и под формулу с учетом релятивисских эффектов.
    
    chersanya
    23.04.2017 22:55
    #10187760
    Здесь же были взяты все формулы в классическом (не релятивистском) случае — при учёте ОТО всё сложнее, конечно, и скорость света не превысится.
    
    maisvendoo
    23.04.2017 23:07
    #10187772
    при учёте ОТО всё сложнее
    
    Не то что сложнее — ужас как намного сложнее. ОТО вообще одна из самых забористых теорий, после теории суперструн)
    
    Однако есть вполне четкие представления о том, как будет вести себя тело, падающее в черную дыру, что описано в популярной форме у Кауфмана в книге «Космические рубежи теории относительности».
    
    Для удаленного наблюдателя тело, падающее в черную дыру будет разгонятся до середины пути, а затем замедляться, пока не замрет в вечной неподвижности на расстоянии
    
    $r = \frac{2G\,M}{c^2}$
    от центра дыры. Это расстояние называется радиусом Шварцильда и его сходство с нерелятивистской формулой 2-й космической скорости лишь забавное совпадение.
    
    В системе отсчета тела падающего в дыру движение будет происходить и после пересечения сферы Шварцильда со скоростью не превышающей скорость света. Как-то так
    
    progman_rus
    24.04.2017 09:36
    #10188080
    еще раз спасибо. Очень интересно.
    Кауфмана читал… давно правда.
    
    progman_rus
    23.04.2017 19:50
    #10187562
    правка:
    
    при бесконечном r ваша П будет стремиться к нулю в пределе.
    
    maisvendoo
    23.04.2017 19:54
    #10187568
    при бесконечном r ваша П будет стремиться к нулю в пределе.
    
    Так и должно быть. Удобным способом представления потенциальной энергии является выбор нулевого её уровня на бесконечности. Все положения ближе характеризуются отрицательной потенциальной энергией.
    
    И вообще, за нулевой уровень потенциальной энергии можно выбрать любую эквипотенциальную поверхность в пространстве где данное поле рассматривается
    
    maisvendoo
    23.04.2017 19:35
    #10187546
    Вот вам и максимальная скорость, при падении из максимально удаленной от планеты
    
    kasyachitche, посмотрите внимательно на Вашу скорость, у Вас ошибка, Вы потеряли двойку и квадратный корень
    
    maisvendoo
    23.04.2017 19:50
    #10187566
    Если интегрировать, то нужно делать это правильно
    
    $\ddot z = -\frac{G\,M}{(R + z)^2}$
    $\frac{dv_z}{dt} = -\frac{G\,M}{(R + z)^2}$
    Учитывая, что $\frac{dz}{dt} = v_z$
    
    $\frac{dz}{dt} \, \frac{dv_z}{dz} = -\frac{G\,M}{(R + z)^2}$
    $v_z \, \frac{dv_z}{dz} = -\frac{G\,M}{(R + z)^2}$
    $\int\limits_0^v v_z\,dv_z = -G\,M\int\limits_{\infty}^0 \frac{dz}{(R + z)^2}$
    $\frac{v^2}{2} = \left\frac{G\,M}{R+z} \right|_{\infty}^0 = \frac{G\,M}{R}$
    или, окончательно
    
    $v = \sqrt{\frac{2\,G\,M}{R}}$
    
    kasyachitche
    23.04.2017 20:23
    #10187604
    Ваша правда. Я лихо положил dt = dz и интегрировал ускорение по расстоянию, полагая, что получу скорость.

kasyachitche
22.04.2017 21:47
#10186524
Почему демонстрация мощной системы ведется на примере решения ДУ? Простейший решатель таких уравнений можно написать самому за несколько минут.

Можно ли попросить развернуть описание, а еще лучше — сравнение с MatLAB? Как там дела с программированием, какие парадигмы поддерживаются, в какой степени? В комментариях писали, что в Scilab нет анонимных функций, так может есть еще какие-нибудь инструменты, облегчающие написание больших программ? Что там с графиками? Насколько широк функционал работы с ними? Какой подход к этой работе?
Насколько удобна и быстра работа с файлами?

А то написали что инструмент мощный, но вот вам пример использования его как калькулятора.
1. sshikov
  22.04.2017 21:53
  #10186534
  Простейший решатель таких уравнений можно написать самому за несколько минут.
  Кстати да. Системы уравнений подобного вида у нас решали все студенты до единого, численно, на Алголе-60, примерно в 1975 году. И по объему кода получалось не сильно больше того, что тут приведено. И решалось это на машине с 16 килословами памяти. Неужто не на чем как-то продемонстрировать прогресс, достигнутый с 1975 года? )))
  1. shbma
    22.04.2017 22:17
    #10186558
    Ну надо же с чего-то начинать). Я взял простейший пример, чтобы показать порядок работы в максимально понятной форме.
    Раз тема интересна сообществу, в ближайшее время сделаю более развернутую статью, где постараюсь ответить на поставленные вопросы по Scilab-у и показать что-нибудь более прогрессивное.
    
    sshikov
    22.04.2017 23:01
    #10186588
    Тема хорошая. С octave бы еще сравнить, хотя бы немножко.
    
    KiloLeo
    23.04.2017 10:13
    #10186920
    +1
    И с Python+Scipy+Matplotlib+…
    Думается, что универсальный питон c библиотеками ни в чём сейчас не уступает специализированным пакетам. Но за счёт универсальности имеет существенное преимущество
    
    sshikov
    23.04.2017 10:57
    #10186968
    Собственно, это всегда так было. Мне сразу вспоминается reduce — универсальный lisp с библиотеками символьной математики.
    
    kasyachitche
    23.04.2017 11:46
    #10187034
    А Simulink?
    У Python есть что-нибудь похожее?
1. maisvendoo
  22.04.2017 23:34
  #10186634
  Почему демонстрация мощной системы ведется на примере решения ДУ? Простейший решатель таких уравнений можно написать самому за несколько минут.
  
  А то написали что инструмент мощный, но вот вам пример использования его как калькулятора.
  
  Что называется «ткнули пальцем в небо».
  
  Теория диффур — один из сложнейших разделов математики. Достаточно навскидку гуглануть, чтобы понять объем статей и монографий, посвященных именно численным методам решения ОДУ, особенно жестких ОДУ, какие например встречаются в динамике поезда.
  
  Так что простейший решатель может и пишется за пять минут, а вот задача численного решения ОДУ и инструментарий, используемый для этого мне лично калькуляторами не кажутся
  1. kasyachitche
    23.04.2017 11:24
    #10187004
    Решить приведенное уравнение — калькуляторная задача. К тому же я сомневаюсь, что функция ode() сложна в своей реализации и подходит для решения сколь-нибудь сложных уравнений или систем. Но это — опираясь на реализацию подобных функций в MatLAB'е.

sshikov
22.04.2017 21:50
#10186530
Формула для сопротивления у вас какая-то сомнительная. Дело в том, что Cx, он обычно экспериментально измеряется, по крайней мере для тел такой сложной формы, как человек, и на таких скоростях (дозвуковых и околозвуковых). Вы, кстати, где его взяли, и какой величины?

А площадь S — это никакая не максимальная, а просто характерная площадь (т.е. грубо говоря, если это Cx для тонкой и плоской пластинки, то это будет площать пластинки, а если цилиндр — то например площать сечения цилиндра.

Ну и кстати — разница в величине g на 30 километрах может и 1%, но если ее влияние проинтегрировать по времени падения — то может получиться совсем не мало.
1. shbma
  22.04.2017 22:39
  #10186572
  Cама-то формула по-моему не подает повода для подозрений, а вот значения упомянутых Вами величин действительно условные.
  
  Cx = 0.35 взял отсюда как для «парашютиста, падающего плашмя, с разведенными в стороны ногами и руками».
  S = 1 взял, что называется, «на глаз»
  
  Если учесть, что герой статьи менял угол наклона тела в процессе полета (вначале падал плашмя, затем под более острым углом), Cx и S по-хорошему должны быть функциями времени полета. Но экспериментальных данных такого рода нет, а придумывать исходя из одного видео затея сомнительная, будет еще хуже.
  1. sshikov
    22.04.2017 23:00
    #10186584
    Ну да, формула норм, а вот циферки да, сильно сомнительны. Впрочем аэродинамика тела человека это штука сложная.
    
    Для S=1 хорошо бы еще понимать, в каких единицах. Квадратный метр?
    
    Cx и S по-хорошему должны быть функциями
    Не совсем так. S это характерная площадь, она не меняется для одного тела. А зависит это все от угла атаки (сопротивление и подъемная сила). График-то несложно получить, и он достаточно простой обычно — но появится еще одна переменная, этот самый угол атаки, которым парашютист может управлять (и еще диф. уравнения).
1. maisvendoo
  22.04.2017 23:15
  #10186606
  +3
  Ну и кстати — разница в величине g на 30 километрах может и 1%, но если ее влияние проинтегрировать по времени падения — то может получиться совсем не мало.
  
  Вы не совсем правы — разница будет несущественна. Выполним её оценку
  
  За начало отсчета примем уровень поверхности Земли. Сопротивление атмосферы пренебрегаем. Ускорение свободного падения зависит от высоты над поверхностью так
  
  $g = \frac{G \, M}{(R + z)^2}$
  Высота прыжка 39 км, радиус Земли — 6378 км, так что можно линеаризовать зависимость ускорения от высоты, разложив вышеприведенную функцию в ряд Тейлора, отбросив малые второго порядка
  
  $g(z) \approx g(0) + g'(0) \, z + ... = \frac{G\,M}{R^2} - \frac{2 \, G \, M}{R^3} \, z + ...$
  Учитывая что
  
  $\frac{G \, M}{R^2} = g_0$
  ускорение силы тяжести на поверхности, получим
  
  $g(z) \approx g_0 \left( 1 - \frac{2 \, z}{R} \right)$
  Составим уравнение движения и проинтегрируем его, считая, что скорость парашютиста v(h) = 0, где h — высота прыжка
  
  $\frac{dv_z}{dt} = - g_0 \, \left(1 - \frac{2\,z}{R} \right)$
  $v_z\,\frac{dv_z}{dz} = - g_0 \, \left(1 - \frac{2\,z}{R} \right)$
  $\int\limits_0^v v_z \, dv_z = -g_0\int\limits_h^0 \left(1 - \frac{2\,z}{R} \right)$
  $v^2 = 2 \, g_0 \, h \, \left(1 - \frac{h}{R} \right)$
  Здесь $2\,g_0\,h = v_0^2$ — квадрат скорости, что получилась бы, считай мы g постоянным. Та что относительное отклонение скорости при этом допущении составит
  
  $\frac{v}{v_0} = \sqrt{1 - \frac{h}{R}} = \sqrt{1 - 0.006115} \approx 0,997$
  Ошибка по скорости составляет 0,3 %. Так что допущение постоянства g вполне справедливо
  1. sshikov
    23.04.2017 09:01
    #10186858
    Да, вы правы. А я как-то подзабыл, что постоянный множитель при интегрировании просто выносится наружу.
    
    Можно было еще проще это показать — просто взять реальную величину g на высоте прыжка, и посчитать все с ней. Чисто на глаз должно получиться отличие между 1% и нулем, что вы и подтвердили.

worldmind
22.04.2017 23:09
#10186596
Говоря о свободном математическом софте думаю стоит упоминать, что не только scilab есть, вот например список.
1. potan
  24.04.2017 13:19
  #10188460
  Туда еще стоит добавить Julia.

eduard93
22.04.2017 23:32
#10186630
Возможно ли получить аналитическое решение с помощью Scilab?
1. shbma
  23.04.2017 08:52
  #10186850
  В Scilab есть возможности для символьных вычислений, но очень ограниченные, в сравнении с упомянутой ниже Maxima.

mrrouter
23.04.2017 00:45
#10186696
Нормально. У нас тоже была эта тема. Даже отчёт по ней где-то сохранился (там тоже были численные решения дифференциальных уравнений, в том числе в частных производных). Кстати, пакет listings для LaTeX признаёт язык скриптов Scilab, что удобно для красивого оформления кода.
К комментарию выше: для этого больше подойдут системы компьютерной алгебры (Maxima и т.д.).

KiloLeo
23.04.2017 10:21
#10186930
Скорость звука на высоте 30 000 равна 295 м/с. Так что звуковой барьер будет перейдён раньше.

Kazancev
23.04.2017 13:14
#10187150
+2
Пару раз открывал Матлаб на лабах по МКЭ в ВУЗ'е, было лень разбираться и качать гибабайтные дистрибутивы — воткнул у себя на компах Scilab. Т.к. преподавателю были важны результаты расчётов, а не софт, то и начал писать МКЭ'шки в Scilab на своих нетбуке и десктопе, на лабах запускал Scilab (можно запускать без установки от админа, скопировав файлы после установки если это Win), показывал результаты и отчитывался.

На мой взгляд, плюсы: 1) мало весит, ~ 150 МБ; 2) для меня понятная справка на смеси English+ местами русский, более понятная чем запутанные портянки матлаб; 3) сижу на Linux + Scilab — полёт нормальный, рабочий вариант Win+Scilab тоже пашет на меня; 4) в Симулинк не работал и не открывал ни разу, работаю в XCos, считаю полёт нормальный, всё моделится ( тут пример с файлом); 5) есть навески-атомы в онлайн-каталоге Atoms, из них пользую вейвлеты, ЦОС, всё остальное для меня есть в Scilab; 6) разрежённые матрицы — люблю, ценю, благодарю разрабов, мегавещь.

P.S. Диплом писал в Scilab (черновик). Зимой парой-тройкой строк написал для себя связку Scilab+CalculiX для параметрического формирования конечно-элементной сетки требуемого объёкта, все операции с ФС работают без напряга с моей стороны, чтение-запись данных идёт, машинка крутится… Работа с кодом выглядит примерно так

.

Инструменты отладки вроде как изменились к лучшему в v6.0.0, но пока не юзал. Симуляция обратного маятника выглядит примерно так:

P.P.S. Для интересующихся есть примеры, в главном окне выбрать Справка-> Примеры. Там можете сравнить функционал матлаба.
1. Kazancev
  23.04.2017 13:21
  #10187160
  Ммм, гибабайт, изобрёл новую величину случайно, звиняйте=)
  1. mrrouter
    25.04.2017 23:13
    #10191152
    Вот и думай: «гигабайт» или «гибибайт» имелся в виду? Хотя куда хуже выглядит употребление слова «диплом» для обозначения дипломной работы и ненужный дефис в слове «конечноэлементный».