Раздел: Документация
0 ... 114 115 116 117 118 119 120 ... 365 % для вектора значений X. Возвращается матрица F, каждый столбец % которой F(:, j) есть f (х (j), у(х(j))). L = length(X) ; % определение количества значений независимого аргумента F = zeros(2, L); % создание нулевой матрицы F подходящего размера % вычисление компонент вектор-функции правой части для всех значений % независимой переменной F(l, :)=Y(2, :); % вычисление первой компоненты F(2, :) = -sin(X); % вычисление второй компоненты Для переключения солвера bvplc па использование векторизованной функции служит свойство vectorized управляющей структуры, которое по умолчанию равно "off. В листинге 6.31 следует заменить подфункцию rside на rsidevect и изменить обращение к солверу: options = bypset(Vectorized, on); sol = bvp4c(Srsidevect, ©bound, initsol, options); Суммарная информация о вычислениях— количество узлов расчетной сетки, максимум невязки и число обращений к функциям правой части системы и граничных условий — выводится в командное окно, если значением свойства stats является on. Класс задач, решаемых солвером bvp4c, включает в себя граничные задачи для систем обыкновенных дифференциальных уравнений с неизвестными параметрами и систем с вырождающимися коэффициентами. В следующих двух разделах демонстрируется решение этих задач на простых примерах. Граничные задачи с неизвестными параметрами Солвер bvp4c так же, как и солверы для задачи Коши и уравнений с запаздыванием, позволяет решать граничные задачи для систем обыкновенных дифференциальных уравнений с параметрами. Однако в случае граничных задач параметры могут быть как известные, так и подлежащие определению в процессе вычислений. Рассмотрим общую постановку граничной задачи с подлежащими определению параметрами. Требуется найти вектор-функцию v(jc) и вектор параметров р, удовлетворяющие на отрезке [а, Ь\ системе дифференциальных уравнений у = /(а, у, р) и граничным условиям: s{y{a)> у(Ь)> р) = ®- Решение таких задач предполагает указание началь- ного приближения как для искомой вектор-функции, так и для всех параметров. Для этого следует инициализировать структуру с начальным приближением, задав три входных аргумента: meshinit — массив узлов начальной сетки, yinit — начальное приближение для искомой вектор- функции (постоянный вектор или указатель на функцию) и parinit — вектор стартовых значений параметров. solinit = bvpinit(meshinit, yinit, parinit) Функция bvpinit создает структуру solinit, поле parameters которой CO-держит начальные приближения для неизвестных параметров. При программировании функций для вычисления правой части системы и граничных условий следует предусмотреть возможность передачи значений параметров в третьем входном аргументе— векторе. Заголовки функций должны иметь вид: function f = rsidepar(x, у, par) И function f = boundpar(ya, yb, par) В случае успешного завершения вычислений солвер bvp4c возвращает структуру sol » sol = bvp4c(@rsidepar, Sboundpar, solinit) Узлы сетки записаны в sol.x, приближение к вектор-функции в sol.у, а найденное значение параметров в sol.parameters. Разберем поиск неизвестных параметров граничных задач на модельном примере с одним параметром. Требуется найти функцию и(х) и число А, удовлетворяющие на отрезке [0, тг] дифференциальному уравнению и" + 2Хи + 2Х2и = 0 и граничным условиям ы(0) = 0, г/(я) = 0. Такого класса задачи часто возникают при определении собственных чисел и собственных функций дифференциальных операторов. Приведенная выше задача, очевидно, имеет тривиальное решение н(х) = 0, не представляющее практического интереса. Как правило, имеется целый спектр параметров, при которых существует ненулевое решение н(лс), определенное с точностью до мультипликативной константы. Наша задача не полностью определена, поскольку если ы(х) — решение, то с-и(л) также является решением для любой константы с. Следовательно, необходимо подчинить искомую вектор-функцию некоторому нормировочному условию, например: и(0) = 1. Введем обозначения у, (л) = н(л), у2 (л) = и{х) и перейдем от дифференциального уравнения второго порядка к системе обыкновенных дифференциальных уравнений первого порядка с соответствующими граничными условиями д(0) = 0. л(я) = 0, Л(0)-1 = 0.(6.12) Примечание В такой постановке эту задачу можно трактовать, как задачу Коши с дополнительным условием на другом конце для определения параметра. При этом нельзя гарантировать существование решения. Создадим теперь функции rsidepar для вычисления правой части системы и boundpar для граничных условий в соответствии с приведенным выше правилом (листинг 6.32). Важно, что хотя граничные условия не зависят от параметра, все равно третий аргумент par должен формально присутствовать в списке входных аргументов функции boundpar. Осталось запрограммировать функцию для определения начального приближения, которое существенным образом влияет на результат вычислений для задач с неизвестными параметрами. В нашем примере хорошим начальным приближением может оказаться У(л) = 5Й1л- (и, соответственно y2(jc) = cos.x), поскольку оно точно удовлетворяет всем трем граничным условиям. Текст функции yinit приведен в листинге 6.33. Теперь сгенерируйте структуру soiinit с информацией о начальном приближении, выбрав сетку с 10 равноотстоящими узлами на промежутке [0, я] и взяв X = 0.8 в качестве стартового значения параметра. Вызовите солвер bvp4c, отобразите графически первую компоненту вектор-функции решения (т. е. искомую функцию «(*)) и найдите параметр Л. Текст файл-функции boundparproblem, осуществляющей перечисленные действия, приведен в листинге 6.33. В качестве выходного аргумента она возвращает значение искомого параметра. Листинг 6.33. Файл-функция для решения граничной задачи j с параметром (6.12) function par = boundparproblem Meshinit = linspace(0, pi, 10); % генерация начальной сетки Pinit =0.8; % стартовое приближение для параметра i генерация структуры с информацией о начальном приближении soiinit = bvpinit(meshinit, ©yinit, pinit); % вызов солвера sol = bvp4c(Srsidepar, ©boundpar, soiinit); % занесение в par найденного значения параметра par = sol. parameters l За* 130 0 ... 114 115 116 117 118 119 120 ... 365
|
