Раздел: Документация
0 ... 113 114 115 116 117 118 119 ... 365 initsol = bvpinit(meshinit, yinit); % вызов солвера sol = bvp4c(grside, @bound, initsol); % использование полей x и у структуры sol для построения решения Щ в sol.х содержатся координаты сетки % в sol.у матрица sol.y(l, :) соответствует значениям функции yl в sol.x(:) Ъ sol.у(2, :) соответствует значениям функции у2 в sol.x(:) plot(sol.x, sol.у(1, :), k.) % вывод графика точного решения для сравнения hold on fplot(@sin, ЕО ll*pi/2]) title{Решение граничной задачи солвером bvp4c) legend(приближенное решение, точное решение) % функция rside для правых частей системы уравнений (6.11) function f = rside(х, у) f = [у(2); -sin(x)]; % функция bound граничных условий function f = bound(уа, yb) i = [уа(1); yb(2) + yb(l) + 1]
Рис. 6.19. Сравнение приближенного и точно! о решений i раничной задачи (6.11) Получающиеся графики приближенного и точного решений исходной задачи приведены на рис. 6.19, они свидетельствуют о том, что для простых задач солвер bvp4c дает хорошие результаты. По умолчанию решение найдено с относительной точностью 1(и и абсолютной 1(н по невязке. Задание точности вычислений, ряда других опций солвера bvp4c и решение более сложных задач обсуждается в следующих разделах. Возможности солвера bvp4c, управление вычислениями Солвер bvp4c позволяет решать граничные задачи для систем обыкновенных дифференциальных уравнений произвольного порядка п у= fix, у), x<=[a,b], g{y(a),y(b)) = 0. Здесь у — неизвестная вектор-функция; g — вектор-функция граничных условий. Правило программирования функции для вычисления правой части системы и вектор-функции граничных условий остается тем же, что и в случае системы второго порядка, рассмотренной в предыдущем разделе. Точность решения контролируется при помощи двух опций управляющей структуры RelToi — для относительной точности и AbsTol — для абсолютной. Сама управляющая структура генерируется функцией bvpset Options = bvpsetСRelToi, 1.0е-05, AbsTol1.Ое-07) и задается в качестве четвертого входного аргумента солвера bvp4c. Абсолютная точность участвует в оценке каждой компоненты вектор-функции невязки г(х) = у(х)-/(дс,5(а:)), где у(х)— приближенное решение, и может быть числом, либо вектором из п элементов. Параметр RelToi используется для оценки интегральной нормы компонент невязки г; на каждом из участков расчетной сетки. Вычисления останавливаются при выполнении условия: \\г( / тах{, AbsTol(i)/RelTol) < RelToi. Возможна ситуация, когда решение граничной задачи найдено, но по каким-то причинам требуется вычислить его с большей точностью. В этом случае в качестве начального приближения для солвера bvp4c допускается указание полученной структуры с приближенным решением. Например, в примере предыдущего раздела было получено решение с установленной по умолчанию точностью: sol = bvp4c(@rside, Sbound, initsol); function F = rsidevect(X, Y) % Вычисление правой части системы дифференциальных уравнений Для нахождения решения, скажем, с относительной точностью Ю-7 достаточно сформировать управляющую структуру options = bvpset(1RelToi1, le-7); и обратиться к солверу sol = bvp4c(@rside, @bound, sol, options); Солвер может прервать вычисления, не достигнув заданной точности, поскольку по умолчанию число узлов расчетной сетки не должно превышать целой части от 1000/«. При появлении такого сообщения следует увеличить максимально допустимое число узлов, установив подходящее значение для параметра Nmax управляющей структуры. В алгоритме солвера bvp4c частные производные функций g и / аппроксимируются конечными разностями. Использование аналитических выражений для них может повысить эффективность вычислений. В случае системы линейных уравнений матрица Якоби правой части системы постоянна и содержащий ее массив задается в качестве значения свойства FJacobian. Для нелинейных систем потребуется 5апрограммировать функцию, так же, как и в случае решения задачи Коши (решение задачи Коши для систем обыкновенных дифференциальных уравнений описано выше в этой главе). Функция, вычисляющая производные вектор-функции g{y„, У/,) граничных условий dg/dya , dg/dyb , должна иметь заголовок function [dgdya, dgdyb] = boundjас(уа, yb) и по заданным у (а) и у(Ь) возвращать векторы dg/dya, dg/dyh. Солвер bvp4c будет вызывать bound j ас для нахождения частных производных, если свойству BCJacobian управляющей структуры установить значение eboundjac. Ускорение работы солвера возможно за счет специального способа программирования функции для правой части системы. В примере предыдущего раздела (листинг 6.31) функция rside по заданной паре значений л и >(<*) вычисляла /(-*, .у(-*)) • Альтернативный вариант интерфейса предполагает возможность вычисления матрицы F = f[X, у(Х)) для заданного вектора значений независимой А переменной и матрицы у(X). Модифицированная функция приведена в листинге 6.32. 1 Листинг 6.32. Векторизованная файл-функция для правой части системы (6.11) \ 0 ... 113 114 115 116 117 118 119 ... 365
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
