Файл: Элементы математического моделирования в программных средах MATLAB 5 и Scilab (Андриевский Фрадков).pdf

введенных обозначений напишем программу вычисления неизвестных t',tm)vm :

S=solve(' v ^ / t ^ l s b / m ' , 'v _ m/(T - t _ l)=a/m' ....

' у . Т Ф ч л / г 1 , Ч _ 1 ' , ' Т ' , ' v . m ' ) t_l=simple(S.t_l)

T=simple(S.T) v_m=simple(S.v.m)

Воператоре solve записана исходная система уравнений

иуказано, относительно каких переменных она должна быть разрешена. Решения системы содержатся в структуре S, поля которой имеют значения класса sym (символьный). Последующие операторы программы осуществляют доступ к полученным результатам с одновременным их упрощением. В результате выполнения программы получим выражения:

[ -y_T*m*(b+a)*2~(l/2)/(m*(b+a)*y_T*a*b)~(1/2)] v_m =

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

1 Отрицательные решения можно трактовать как движение в противоположном направлении и обратном времени.

40

а вызов latex (tl) приводит к строке символов, которая после обработки текстовым редактором IATjtX[57] дает выраже-

Теперь перейдем к решению уравнения (1.4) при заданных начальных условиях и найденном через (1.13) управлении. Начнем с символьных вычислений. Для этого продолжим программу следующими операторами:

vy=dsolve('m*Dv=b','Dy=v','у(0)=0','у(0)=0'); v=simple(vy.v)

y=simple(vy.у)

v _ t l = s u b s ( v , ' t ' , t _ l ( l ) ) y _ t l = s u b s ( y , ' t ' , t _ l ( l ) )

В этом фрагменте первая строчка содержит обращение к процедуре символьного (аналитического) решения уравнения (1.4), приведенного предварительно к системе уравнений первого порядка относительно переменных v(t) и s(<), имеющей вид

Символом D в программе обозначается операция дифференцирования. Переменные v, у содержат найденные при = 6 решения. Оператором subs выполняется подстановка в них значения t = t', найденного выше. В формате IATj?X[57] получаются выражения

v

Продолжим вычисления. Рассмотрим следующий интервал времени t 6 [t'}T]. На нем u(t) = —а. Проинтегрируем нашу систему на этом участке. Для этого выполним операторы:

vy2=dsolve(,m*Dv=-a', 'Dy=v', 'v(tl)=vl' , 'y(tl)=yl') v2=simple(vy2.v)

y2=simple(vy2.y)

vT=simple(subs(v2,{'vl','t1','t'>,{v _ tl,t . l(1),T(1)>))

Оператором dsolve находится решение (1.4), для которого в момент t = выполнено v(ti) = vX) y(t 1) = Значения tX)Vi, yx подставляются оператором subs по результатам вычислений на предыдущем участке. Получим

41

v2 =vl - a*(t - tl)/m

y2 =yl+vl*t-vl*tl-l/2*a*(t~2+tl~2-2*t*tl)/m vT =0

Равенство нулю значения скорости на правом конце интервала соответствует условию поставленной задачи.

2.2.2.Применение процедуры численного интегрирования дифференциальных уравнений

Обратимся теперь к численному решению уравнения (1.4). Сначала приведем программу моделирования для более общей ситуации, а именно не будем учитывать, что рассматривается линейное уравнение с постоянными параметрами. Воспользуемся одной из процедур интегрирования обыкновенных дифференциальных уравнений, входящих в систему MATLAB, а именно процедурой ode45, в которой используется метод Рунге-Кутта 4-5-го порядков [74].

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

Головная программа

clear close a l l

global a b m t l a=2e3;

b=5e3;

1 В пакете MATLAB имеется несколько процедур численного интегрирования дифференциальных уравнений. Для каждой из них может быть установлено значительное число параметров (опций). Более подробную информацию можно получить из документации, или выполнив команду help ode45. В нашем примере используются значения параметров, принятые по умолчанию.

42

m=le3;

у.Т=10;

tl=l/(b+a)*2~(l/2)#(m*(b+a)*y_T#a*b)л (1/2)/Ь T=y_T*m*(b+a)*2~(1/2)/(m*(b+a)*y_T*a*b)~(1/2)

global a b m t l dx=zeros(size(x)); if t<tl

u=b; else

u=-a; end

dx(l)=u/m;

dx(2)=x(l);

Рассмотрим эти программы подробнее. Прежде всего заметим, что если имя rn-файла головной программы безразлично (более того, ее операторы могут быть введены непо-

[35, 74, 84] и помещается в рабочем директории (папке) под определенным именем. В рассматриваемом примере текст этой программы содержится в файле fpl_2.m. Имя программы (совпадающее с именем файла), записано в ее заголовке

function dx=fpl . 2(t,х);

и, что более существенно - указывается в головной программе при обращении к процедуре интегрирования:

Ct,x]=ode45( , fpl . 2 , ,[0, Т],х0);

43

Головная программа начинается операторами clear и close. Эти операторы необязательны, но с них рекомендуется начинать головную программу. Первый из них освобождает рабочую область от всех переменных, которые могут там оказаться после предыдущих вычислений.1 Оператор close закрывает графические окна, если они были открыты перед выполнением программы. Без этого каждый "запуск" программы будет сопровождаться увеличением их числа. В последующих примерах операторы clear и close указываться не будут. В операторе global перечисляются глобальные переменные, значения которых доступны как в головной, так и в вызываемой программах без передачи через параметры процедуры. Такой же оператор имеется и в подпрограмме fpl_2. В следующих операторах (а=2еЗ; Ъ=5еЗ; т=1еЗ; у_Т=10;)

(в этом примере вектор переменных состояния х = [v,y] ). Далее производится обращение к процедуре ode45, в котором указаны входные параметры: имя процедуры-функции вычисления правых частей уравнений; вектор [0,Т], содержащий начальный и конечный моменты времени, для которых нужно проинтегрировать уравнение; начальное значение состояния системы. Выходными величинами являются массив моментов времени t, для которых получено решение, и массив результатов интегрирования х. Массив t является вектор-столбцом, число элементов nt которого определяется автоматически процедурой интегрирования. Значение nt можно получить, выполняя по окончании вычислений опе-

раторы l e n g t h ( t ) или s i z e ( t ) . В данном примере nt = 73. Массив результатов имеет размер ntxn, где п - порядок ин-

тегрируемой системы (в нашем случае п = 2). Каждая строка массива х относится к соответствующему моменту вре-

1 Естественно, если программист хочет в дальнейшем использовать эти переменные, то оператор clear выполнять нельзя.

2 Заметим, что в данном примере выражения для этих переменных получены копированием результатов символьных вычислений, описанных выше.

44

смотрено ниже). Далее в головной программе выполняется вывод результатов на график оператором p l o t ( t , x ) . В той же строке через запятую размещен оператор grid, задающий

да заголовка графика и надписи по оси аргументов (х-ось). Оператор legend( ' v ' , 'у' ,0) служит для вывода легенды - связи имен переменных с линиями на графике. Цифра "О" в последнем аргументе оператора legend означает, что для меток должно автоматически выбираться "наиболее удобное" (свободное от графика) место. 2

Результат работы программы показан на рис. 2.3, а. Время интегрирования на компьютере с процессором Pentium 166 ММХ составляет 0.15-г 0.17 с. Рассмотрим теперь программу

v. м/с; у. и

Рис. 2.3. Результаты интегрирования и ошибки вычислений.

fpl_2 (стр. 43). Состав формальных параметров этой подпрограммы фиксирован и определяется требованиями процедуры ode45: это входные значения момента времени t и переменных состояния х и выходной массив вектора правых ча-

1Перечень свойств объекта и их значений можно получить оператором set (в нашем случае - вызовом set(gca)) [35, 74, 84]. Отметим, что системе MATLAB 5.3 имеется удобный интерактивный редактор графиков, использование которого позволяет упростить текст программы.

2При выводе на терминал графики имеют разные цвета. Это отличие пропадает при черно-белой печати.

45