Раздел: Документация
0 ... 100 101 102 103 104 105 106 ... 365 Решите задачу, используя, например, солвер odeii3. Входными аргументами солверов в самом простом случае являются: указатель на функцию (или ее имя в апострофах), вектор с начальным и конечным значением времени наблюдения за процессом и вектор начальных условий. Выходных аргументов два: вектор, содержащий значения времени, и матрица значений искомых функций в соответствующие моменты времени. Значения функций расположены по столбцам матрицы, в первом столбце— значения первой функции, во втором— второй ит. д. В силу проделанных замен ух — у, Уг~У-> первы& столбец матрицы содержит как раз значения неизвестной функции у(/), входящей в исходное дифференциальное уравнение, а второй столбец — значения ее производной. Как правило, размеры матрицы и вектора достаточно велики, поэтому лучше сразу отобразить результат на графике. Применение солвера для нахождения решения при ( < 15 и визуализация результата продемонстрированы на примере файл-функции solvdem, приведенной в листинге 6.15, где функция oscil реализована как подфункция. Листинг 6.15- Файл-функция soivdan для решения дифференциального уравнения (6.4) function solvdem % формирование вектора начальных условий УО = [1; 0]; % вызов солвера от функции oscil, начального и конечного % момента времени и вектора начальных условий [Т, У] = odell3(@oscil, [0 15], У0); % вывод графика решения исходного дифференциального уравнения % (маркеры - точки, линия - сплошная) plot(T, Y(:, 1), г.-) % вывод графика производной от решения исходного % дифференциального уравнения (маркеры - точки, линия - пунктир) hold on plot(T, Y(:, 2), *:..-) % вывод пояснений на график title(Решение (\ity) \prime\prime+2{\ity} \prime+10{\ity) =sin(\itt)) xlabeMUtt) ylabel { {Uty} , {\ity} \prime 1) legend(1 координата, скорость, 4) grid on hold off % подфункция вычисления правых частей уравнений function F = oscil(t, у) F = [у(2>; -2*у(2) - 10*у(1) + sin(t)]; В результате выполнения файл-функции solvdem на экран выводятся графики, изображенные на рис. 6.8, которые отражают поведение координаты точки и ее скорости в зависимости от времени. Из графика видно, что приближенное решение и его производная удовлетворяют начальным условиям, колебание происходит в установившемся режиме, начиная с t - 5. Примечание Солверы могут найти приближенное решение для заданных значений независимой переменной. Для этого следует в качестве второго входного аргумента указать вектор с этими значениями, упорядоченными по возрастанию или убыванию. Решение у "+2у 410y=sirtf  -1.5 -2 -2.5 - координата --- скорость 1D IS Рис. 6.8. Решение дифференциального уравнения (6.4) Заметим, что интерфейс солверов допускает обращение к ним с одним выходным аргументом: » sol = odell3(@oscil, [0 15], Y0); Такой вызов солвера приводит к образованию структуры sol с информацией о полученном приближенном решении. Поле х структуры sol содержит вектор-строку значений независимой переменной, а поле у— матрицу из соответствующих значений компонент решения, записанных в строки (работа со структурами данных подробно описана в главе 8). Проще говоря, soi.x эквивалентно Т, a sol.у— y1 для рассмотренного выше обращения к солверу с двумя выходными аргументами т и y в файл-функции solvdem, текст которой приведен в листинге 6.15 (напомним, что апостроф означает транспонирование). Если бы в solvdem солвер был вызван с единственным выходным аргументом, то для построения графиков компонент приближенного решения следовало бы использовать команды (проверьте!): plot(sol.х, sol.yd, :), г.-) plotfsol.x, sol.y(2, :), k.:) Мы привели альтернативный вариант интерфейса солверов с одним выходным аргументом, поскольку он позволяет ответить на вопрос: как получить значение решения в любой точке отрезка. Действительно, солверы вычисляют приближенное решение задачи Коши в конечном числе точек на заданном отрезке, причем координаты точек выбираются солвером или задаются пользователем. Для определения значений компонент решения в промежуточных точках придется прибегнуть к интерполяции. Эффективную интерполяцию осуществляет функция deval, использование которой требует структуры sol с информацией о решении. Самый простой способ вызова функции deval >> sxint = deval(sol, xint) предполагает указание вектора xint с координатами точек, в которых следует вычислить решение. Найденные значения компонент построчно записываются в выходном аргументе, т. е. в sxint (:, i) хранятся значения всех компонент искомой вектор-функции в точке с координатой xint(i). Например, для вычисления компонент решения на отрезке [0, 5] в равноотстоящих точках с шагом 0.02 и построения графиков достаточно выполнить команды: » xint = 0:0.02:5; » sxint = deval(sol, xint); » plot(xint, sxintd, :), xint, sxint(2, :)) Интерполяция возможна только для некоторых компонент решения, в этом случае необходимо указать их номера в векторе idx и включить его в список входных аргументов » sxint = deval(sol, xint, idx) 0 ... 100 101 102 103 104 105 106 ... 365
|
