Математическое моделирование: аналитические и вычислительные методы Текст научной статьи по специальности «Математика»

Кондратьев Александр Сергеевич, Ляпцев Александр Викторович

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ: АНАЛИТИЧЕСКИЕ И ВЫЧИСЛИТЕЛЬНЫЕ МЕТОДЫ

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

Следует отметить, что в последнее время такая точка зрения теряет своих поклонников, и все более и более ясным становится положение, что достаточно сложные модели приводят к уравнениям, не решаемым в аналитическом виде [1]. Наблюдающееся в последнее время бурное развитие численных методов, в частности, вычислительной физики, является наглядным отражением этого обстоятельства.

Развитие численных методов рассмотрения математических моделей должно сопровождаться развитием качественных методов исследования уравнений, что позволяет установить общую картину явлений без

решения соответствующих уравнений и способы контроля правильности проводимого численного решения. Как отмечает Р.Фей-нман, «Грядущая эра пробуждения человеческого разума приведет к пониманию качественного аспекта уравнений. Сейчас мы еще не способны на это» [2].

Обучение основам математического моделирования должно отражать обе указанные тенденции - развитие как умений качественного анализа, так и умения нахождения численных решений. И здесь очень важно с самого начала установить правильный акцент на то, что не всегда аналитическое решение является предпочтительным: иногда оно оказывается неоправданно громоздким и трудоемким по сравнению с численным методом решения, при этом не приводя к какой-то новой полезной информации. В качестве примера можно указать на «классическую» задачу о соскальзывании шайбы с поверхности полусферы. В школьном варианте этой задачи, когда трение отсутствует, имеется прекрасная возможность качественного исследования характера движения шайбы с доказательством отрыва ее от поверхности полусферы [3].

А °

^ кагес&$е примера, мафНа цкауг&ь На «&лассигесЩ(а» ¿а^агу а саскалъ^м&аНии майби па паёерхНасЛи пам^срери.

При учете трения движение шайбы удобно рассматривать с помощью численных методов решения получающихся дифференциальных уравнений [4]. Однако при определенной постановке задачи можно получить и аналитическое решение. Проектирование векторного уравнения движения шайбы

mg + N + Fmp = ma

(1)

на касательное и нормальное (радиальное) направление приводит к двум уравнениям:

(2)

v

mg cos в - N = m— , R

du

mg sin в Fmp = m — . (3)

Здесь m - масса шайбы, R - радиус полусферы, g - ускорение свободного падения, в - угол, образуемый с вертикалью радиус-вектором шайбы с началом в центре полусферы. Сила сухого трения скольжения F связана с нормальной силой N ре-

mp ^ f

акции опоры соотношением:

FmP = m,

(4)

где т - коэффициент трения скольжения.

Уравнения (2) и (3) справедливы, пока шайба скользит по поверхности полусферы. Отрыв от поверхности происходит при обращении в нуль силы N : N = 0. Из уравнения (2) следует, что соответствующий угол вг определяется соотношением:

в u

cose1 = —1

1 gR

(5)

где ьг - скорость шайбы в точке отрыва.

Для определения нужно решить систему уравнений (2) и (3) при условии (4). Это можно сделать, если представить производную du/dt в виде:

du _ du q _ du u _ 1 d . 2 ~dt _ ~dd' _ 1в' ~R _ 2Rdß(u ).() Здесь учтено, что пока шайба движется по поверхности полусферы, справедливо соотношение u _ RQ . Подставляя в уравнение (3) соотношение (6) и соотношение (4), в котором N выражено с помощью уравнения (2), приходим к уравнению:

d

Это обыкновенное дифференциальное уравнение вида:

у;- 2ту = /(в), (8)

где у = и2, и /(в) - заданная функция угла в, определяемая правой частью уравнения (7).

Общее решение уравнения (8) складывается из общего решения однородного уравнения

ув- 2ту = 0 (9)

и какого-либо частного решения неоднородного уравнения (8). Решение уравнения (9) дается выражением:

у(в) = Сехр(2тв), (10)

а решение неоднородного уравнения может быть получено, например, методом вариации произвольной постоянной С в (10). В результате получаем для общего решения уравнения (7) выражение:

28Я

V2 (в)

((2m -1) cos в —

de

(V2) - 2ßV2 = 2gR(sin в - ß cos в). (7)

1 + 4т

- 3msinв) + С ехр(2цв). (11) где С1 - произвольная постоянная. Подчеркнем еще раз, что уравнения (2) и (3) справедливы только для N > 0 (когда шайба находится на поверхности полусферы).

Для определения С1 необходимо задать начальное условие. Рассмотрим случай, когда шайба начинает скользить без начальной скорости (и = 0) при значении угла в0, определяемого соотношением:

в0 = агС^ т. (12)

Значение в0 соответствует предельному углу, при котором шайба удерживается неподвижной на поверхности полусферы в результате действия сил трения (мы считаем, что максимальное значение силы трения покоя равно силе трения скольжения). Таким образом, начальное условие имеет вид:

и(в0) = 0.

Реализацию подобного процесса опытным путем можно осуществить, например, медленно подталкивая шайбу, находящуюся на вершине полусферы, вниз до начала скольжения. Численный расчет при других начальных условиях рассмотрен в [4]. При рассматриваемых начальных условиях, так

Кондратьев А.С., Ляпцев А.В.

же как и в случае отсутствия трения, легко показать, не решая уравнений, что шайба обязательно оторвется от поверхности полусферы.

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

с =--2R2 ((2m2 - 1)cose0 -

i + 4m

- 3m sine0)exp(-2me0). (13)

Выражение (11) для и2(в) при учете (12) и (13) принимает вид после весьма громоздких, но не сложных преобразований:

U2(в) =-2gR^— (2cose -

4 - 3 cos2 e0

- 3 cos в0 cos(в - в0) +

+ ехр(2(в - во )tgeo )coseo). (14)

Теперь с помощью соотношений (5) и (14) можно получить уравнение для определения угла в1, при котором происходит отрыв шайбы от поверхности полусферы:

- 3 008 в0 С08(в1 - в0) +

+ ехр(2(;-во)18во)со8во). (15)

Возникает естественный вопрос: что проще - непосредственно решать численно систему уравнений (2)-(3) с условиями (4) и (5) или совершать достаточно утомительные и громоздкие аналитические преобразования и в результате прийти к сложному трансцендентному уравнению (15), которое все равно придется решать численными методами? Заметим, что если постановка задачи такова, что необходимо найти зависимость угла в от времени, то необходимо провести интегрирование

в Ив'

' = к ¡Ш ° с(;), (16)

в0

где функция Е(в) есть квадратный корень из выражения в правой части равенства (14), а затем найти функцию, обратную функции О(в). Результат и того, и другого действия не выражается через элементарные функции, то есть аналитические вычисле-

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

Заметим также, что при обучении математическому моделированию целесообразно использовать интегрированные среды, такие как Maple, Mathematica, MATLAB, MatCad, где имеются стандартные и достаточно хорошо работающие процедуры решения систем дифференциальных уравнений, которые могут быть достаточно легко освоены как студентами, так и школьниками старших классов.

Могут ли дать конкретно для данной задачи аналитические методы что-либо полезное? Чтобы ответить на этот вопрос можно попробовать численно решить систему дифференциальных уравнений (2)-(4) при условии и(в0) = 0. В зависимости от реализации численного метода можно получить для скорости либо тождественный нуль, либо что-то не отличающееся от нуля в пределах заданной погрешности вычислений. Чтобы получить нечто, соответствующее реальному явлению, необходимо незначительно изменить условия, например, вместо в0 взять начальный угол в0 = (1 + £)в0, где 0 < £ << 1. Насколько будут различаться расчеты при различных значениях £ и можно ли подойти к пределу при £ ® 0?

Результаты численных расчетов, полученных при решении системы дифференциальных уравнений, приведены на рис. 1, 2. При проведении расчетов сделано масштабное преобразование времени, введена безразмерная переменная t' = t^g / R , что означает измерение времени в единицах ■Jr/ g . Оба рисунка соответствуют одному значению m = 0,1, но разным значениям £ : £ = 0,01 и £ = 0,001. Фазовые диаграммы для обоих случаев, то есть зависимости скорости от угла, определяемые равенством (14), выглядят практически одинаково и приведены на рис. 3. Для сравнения численных методов и аналитических приведена сплошная кривая, соответствующая численному расчету, и точки, соответствующие расчету по формуле (14). Как видно из рис. 1 и 2,

с уменьшением значения е, то есть с приближением начального угла к 0О, время движения шайбы увеличивается, что соответствует интуитивным представлениям о ее движении. Ответить на вопрос, есть ли конечный предел времени движения при е ® 0 численный расчет не может в силу того, что точность расчета всегда ограничена.

Ответить на этот вопрос можно, проанализировав интеграл (16), а именно рассмотрев предел подынтегрального выражения при в' ® в0. Для этого следует разложить правую часть выражения (14) в ряд по степеням малого параметра Дв = в - в0. Из того, что это выражение соответствует начальным условиям и(в0) = 0, следует, что нулевой член разложения равен нулю. Однако несложные расчеты показывают, что и первый член разложения (слагаемое, линейное по Дв) равен нулю, и лишь квадратичный член разложения отличен от нуля:

2gR

и2 (в)

0 -етв0 (1 / 2 + 4 - 3 ^ в0

+ 2 tg2 в0 )(в-в0 )2.

(17)

t = R

■ йв' (в'-вй)

+1

йв'

(19)

Из этого выражения следует, что при значениях в, близких к в0, скорость является линейной функцией от разности

в- в

и(в) » ^0 (в-в0), (18)

где - константа, определяемая равенством (17). Заметим, что эта линейная зависимость «на глаз» видна на рисунке 3, однако очевидно, что строго подобный вывод можно сделать лишь, исходя из аналитического рассмотрения.

Подставляя выражение (18) в интеграл (16), получим, что вблизи в0 можно аппроксимировать выражением:

где значение в1 следует брать настолько близким к значению в0, чтобы можно было пренебречь членами разложения функции и2(в) по Дв более высокими, чем (Дв)2.

Рис. 1

Рис. 2

Рис. 3

в

Кондратьев A.C., Ляпцев A.B.

Второй интеграл в правой части равенства (19) является конечным, а первый - расходящимся. Если же в качестве нижнего предела интегрирования вместо значения в0 взять значение (1 + е)в0, то вычисление первого из интегралов приводит к выражению:

1

In

в -в,

0

ево

Это означает, что при стремлении е к нулю время движения стремится к бесконечности по логарифмическому закону.

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

Литература

1. Самарский A.A., Михайлов А.П. Математическое моделирование. М., 2001.

2. Feinmann R.P. The Pleasure of Finding Things Out. Perseus Publishing, Cambidge, Massachusetts, 1999.

3. Кондратьев A.C., Ситнова E.B. Парадоксальность физического мышления. СПб.: РГПУ им. А.П. Герцена, 2007.

4. Кондратьев A.C., Ляпцев A.B. Физика. Задачи на компьютере. М.: Физматлит, 2007.

Кондратьев Александр Сергеевич, академик РАО, доктор физико-математических наук, профессор кафедры «Методика обучения физике» РГПУ им. А.И. Герцена,

Ляпцев Александр Викторович, доктор физико-математических наук, заведующий кафедрой «Методика обучения физике» РГПУ им. А.И. Герцена.