Решение краевых задач. Лабораторная работа №2б (1253748), страница 3
Текст из файла (страница 3)
В противном случае он становится вычислительно неустойчивым,даже если решение задачи (2) зависит от входных данных «умеренно».При решении уравнений F()=0 методом деления отрезка пополам, задают 0 и 1 так,чтобы разности y(1, ) - Y и y(1, ) - Y имели разные знаки. Далее, полагая = 0 1 , и011122вычисляют y(1,2). Затем вычисляют 3 по одной из формул: 3= 1 2 2 вили 3= 02зависимости от того, имеют ли разности y(1,2) - Y1 и y(1,1) - Y1 соответственно разные или2_____________________________________________________________________________ 7Деменков Н.П.
Решение краевых задачМосковский государственный технический университет им.Н.Э.БауманаКафедра “Системы автоматического управления”__________________________________________________________________________________одинаковые знаки. Далее вычисляется y(1,3). Процесс продолжается до тех пор, пока небудет достигнута требуемая точность | y(1,n) - Y1 |<.В случае использования для решения уравнения F()=0 метода Ньютона задают 0, азатем последующие n вычисляют по рекуррентной формулеn+1= n - F ( n ) , n=0,1,…F ' ( n )Производная F’(n) может быть вычислена по одной из формул численногодифференцирования, например, первого порядка аппроксимации:F’(n) F ( n h) F ( n ) .hПример 1. Пусть требуется решить краевую задачу следующего видаtx'' - x' = 3t2; x(1) = 2; x(2) = 9.Представим данную задачу в виде системы уравнений первого порядка: x1 ' x2 ,x2 x2 ' 3t ,t x1 1 2, x1 2 9.(1)Далее сведем данную задачу к задаче Коши, введя параметр , равный неизвестномузначению x2(1). x1 ' x 2 ,x2 x 2 ' 3t ,t x1 1 2, x 2 1 .Решив полученную систему с фиксированным параметром , мы получим значение x2(2),вообще говоря, отличающееся от истинного.Для корректировки параметра рассчитаем его новое значение по формуле new old x1 2calc x1 2.x2 2calcЗдесь x1(2)calc – полученное в результате расчета значение x1(2).
Затем снова решаемсистему дифференциальных уравнений и т.д.Процесс расчета продолжается до тех пор, пока не будет выполнено условие:new – old , где – заранее заданная точность расчета.Замечание. Для того чтобы найти , при котором выполнено граничное условие в точкеt = 2, удобно ввести еще два уравнения, продифференцировав исходную систему по параметру и введя еще две дополнительные переменные x3 x1x, x4 2 .Для решения краевой задачи создадим в рабочей области MATLAB m-файл boundary.m,,который является функцией переменных t и x. Для создания файла воспользуемся редакторомMATLAB Editor/Debugger, который вызывается из основного меню File – New – M-File.
Текстфайла:function boundary% Задание начальных условийx=zeros(2,1);n=0;ips=10^(-8);alpha=2.5;alpha_old=3.5;_____________________________________________________________________________ 8Деменков Н.П. Решение краевых задачМосковский государственный технический университет им.Н.Э.БауманаКафедра “Системы автоматического управления”__________________________________________________________________________________% Условие окончания процесса интегрированияwhile abs(alpha-alpha_old)>ipsx0 = [2 alpha];t0 = 1;tk = 2;x1end=9;% Интегрирование[t,x]=ode45(@task,[t0 tk],x0);x1 = x(length(x),1)x2 = x(length(x),2)x1endCalc=x1; x2endCalc=x2;alpha_old=alpha;alpha=alpha_old-(x1endCalc-x1end)/x2endCalc;n=n+1;end% Печатьfigure;hold on;plot(t,x(:,1),'-r',t,x(:,2),'-g');grid on;nalpha% Вычисление правых частейfunction dx=task(t,x)dx= [x(2);...3*t+x(2)/t];endendПосле запуска и соответствующего числа итераций программы можно увидетьграфик решения:1211109876543211.11.21.31.41.51.61.71.81.92Для контроля выполнения программы выведены конечные значения x1, x2, n и полученноезначение .ПРИЛОЖЕНИЕ 3Стандартные процедуры для решения краевых задач в MATLABДля решения краевых задач можно использовать специальные программы в MATLABbvp4c или bvp5c.Функция bvp4c решает краевую задачу для системы обыкновенных дифференциальныхуравнений y’=f(y,x).Решение y=y(x), удовлетворяющее на отрезке [a,b] этому уравнению и граничнымусловиям, наложенным на значение функции и/или ее производной на концах отрезка, ищется вформе сеточной функции.
Отрезок делится точками a=x1<x2<…<xn=b (не обязательноравные) и каждой точке xi ставится в соответствие yi. Исходя из начальных значений x1 и y1,путем аппроксимации производной y’=dyв каждой точке сетки разностным отношениемdx_____________________________________________________________________________ 9Деменков Н.П. Решение краевых задачМосковский государственный технический университет им.Н.Э.БауманаКафедра “Системы автоматического управления”__________________________________________________________________________________y, строится система уравнений, из которой находятся значения xi.
В процессе ее решенияxсетка xi может перестраиваться (в частности, сгущаться), при решении используетсяматрица-якобианf.xПростейшая форма обращения к функции bvp4c:sol=bvp4c(odefun,bcfun,solinit).Здесь:odefun- функция, вычисляющая вектор правых частей;bcfun - функция, вычисляющая вектор граничных условий, две компоненты которогопредставляют собой выражения, обращающиеся в 0 в точках a и b соответственно;аргументами функции bcfun являются ya и yb – векторы решения y в точках a и b;solinit – выходная структура функции bvpinit, с помощью которой задаются начальныезначения xi и yi .solinit=bvpinit(xinit,yinit).При этом заполняются два поля: solinit.x=xinit и solinit.y=y(i,:).Здесь xinit – вектор-строка a=xinit(1)<xinit(2)<…<xinit(n)=b, yinit – гипотетическиезначения для y(i), которые могут задаваться в одной из двух форм: в виде вектора, каждаякомпонента которого yinit(i) копируется в качестве гипотетического решения для всех точексетки, т.е.
y(i,;)=yinit(i) или в виде функции в форме y=quess(x), где x – любая точкаотрезка[a,b], y – вектор, длина которого равна порядку системы дифференциальныхуравнений. Для каждой точки сетки x(i) вычисляется вектор гипотетического решенияy(i,:)=quess(x(i));sol – структура, аналогичная solinit, содержит решение краевой задачи; кроме полей sol.x иsol.y она имеет поле sol.yp, в котором содержатся значения производной решения (sol.y)’ вточках sol.x.Полная форма обращения к функции bvp4c:sol=bvp4c(odefun,bcfun,solinit, options,P1,P2,…),где options – аргумент, позволяющий задавать различные управляющие параметры, P1,P2 ит.д.- дополнительные аргументы для вычисления odefun и bcfun. Управляющие параметрызадаются путем обращения к функции bvpset, аргументом которой являетсяпоследовательность пар вида < ‘параметр’, значение >, например,options=bvpset(‘FJacobian’,@FJac).Параметр FJacobian задает в аналитическом виде якобиан [ f ], где f – функция,yвычисляемая в odefun.
Параметр BCJacobian задает в аналитическом виде два якобиана [ bc ]y aи [ bc ], где bc(ya,yb) – функция, вычисляемая в bcfun.ybНеизвестный параметр или вектор неизвестных параметров вводится с помощьюфункции bvpinit, форма обращения к которойsolinit=bvpinit(xinit,yinit, parametrs).Здесь parametrs – гипотетическое значение неизвестного параметра (вектора), solinit –структура, в которой кроме полей solinit.x=xinit и solinit.y=y(i,:) заполняется еще одно полеsolinit.parametrs=parametrs.Пример 2.
Решить задачу (1), используя функцию bvp4c.Составим m-файл, содержащий функцию, вычисляющую правые части:function dxdt=exampl(t,x)dxdt=[x(2);...3*t+x(2)/t];a также m-файл, содержащий функцию, задающую граничные условияfunction res=border(xa,xb)res=[xa(1)-2 xb(1)-9];_____________________________________________________________________________ 10Деменков Н.П.
Решение краевых задачМосковский государственный технический университет им.Н.Э.БауманаКафедра “Системы автоматического управления”__________________________________________________________________________________В качестве начальных значений координаты t выберем с помощью функции linspaceпять равноотстоящих точек на интервале [1,2].tinit=linspace(1,2,5);В качестве начальных приближений вектора решений в этих точках выберемxinit=[0 2]. При таком выборе вектора xinit поле структуры solinit.x будет заполненокопиями исходного вектора. Такой выбор начального приближения может оказаться не оченьудачным.
Альтернативный вариант состоит в задании вспомогательной функции,вычисляющей эти компоненты. Начальное приближение влияет на число итераций. Выборначального вектора или вида вспомогательной функции полностью зависит от интуициипользователя.Сформируем структуру solinit, используемую при решении задачи:solinit=bvpinit(tinit,xinit).Для решения задачи выполним командуsol=bvp4c(@exampl,@border,solinit).Если попытаться построить решение краевой задачи в виде графика с помощьюкомандыplot(sol.t,sol.x(1,:)),результат получится не очень гладкий – слишком мало точек. Эту трудность можнопреодолеть с помощью дополнительной функции deval, которая, используя информацию,содержащуюся в структуре sol, строит интерполяционный сплайн Эрмита, для заданноговектора пробных точек tt, который может содержать сколь угодно точек для обеспечениятребуемой гладкости результата.
Обращение к функции bvpval имеет вид:xx=deval(sol,tt),где xx – значение сплайн-функции в пробных точка. Вектор, содержащий по умолчанию 100пробных точек, создадим командойtt=linspace(1,2)._____________________________________________________________________________ 11Деменков Н.П. Решение краевых задач.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
