Материал: Matlab. Практический подход. Самоучитель

Самоучитель Matlab

и конечное значение для элементов списка и количество этих элементов. Создается список равномерно распределенных элементов. Например, командой linspace(0,1,10) создается список из десяти элементов с первым элементом 0 и последним элементом 1, шаг дискретности равен 0.1. Вся команда по созданию структуры с начальным распределением для искомых функций имеет вид IC=bvpinit(linspace(0,1,10),[-2 0]). Результат записывается в переменную IC, а вторым аргументом функции bvpinit() передается список [-2 0]. В этом случае для первой функции во всех узловых точках используется значение -2, а для второй функции во всех узловых точках используется значение 0. После этого командой sol=bvp4c('MyBVP','MyBC',IC) решаем краевую задачу. Вся последовательность команд выглядит так:

>>IC=bvpinit(linspace(0,1,10),[-2 0]);

>>sol=bvp4c('MyBVP','MyBC',IC);

>>plot(sol.x,sol.y(1,:),'r-','LineWidth',2)

>>grid on

>>hold on

>>plot(sol.x,-3.*sinh(sol.x)./sinh(1)+2.*sol.x,'b:','LineWidth',3)

Результат вызова функции bvp4c() записывается в переменную sol. Это структура, и у нее несколько полей. Поле x (полная ссылка на это поле имеет вид sol.x) содержит набор узловых точек. Поле y (полная ссылка на это поле имеет вид sol.y) содержит значения искомых функций в узловых точках. Значенияпервойфункциисодержатсявстрокеsol.y(1,:)(функция y(x)), для второй функции значения содержатся в строке sol.y(2,:) (функция z(x) ).Поэтомучтобыпостроитьграфикдлявычисленнойфункции,можноиспользовать команду plot(sol.x,sol.y(1,:),'r-','LineWidth',2)

(красная линия толщины 2). Для сравнения строится кривая аналитического решения командой plot(sol.x,-3.*sinh(sol.x)./sinh(1)+2.* sol.x,'b:','LineWidth',3) (синяя пунктирная кривая). На рис. 6.26 показан результат создания кривых для аналитического решения и решения, вычисленного числовыми методами.

Как и во всех предыдущих случаях, совпадение числового и аналитического результатов более чем приемлемое.

Для разнообразия решим еще одну краевую задачу, но уже с использованием функции bvp5c() и с несколько иным способом передачи аргумен-

тов. Задача формулируется так: решить дифференциальное уравнение y′′(x) −y′(x) = 0 с граничными условиями y(0) = −1 и y′(1) −y(1) = 2 .

Задача имеет решение y(x) = exp(x) −2 . Общая схема решения стандартная: сводим задачу к системе уравнений первого порядка, создаем функции для системы уравнений и граничных условий, начальное распределение и затем с помощью встроенной функции находим решение.

266

Глава 6. Интегрирование и дифференциальные уравнения

Рис. 6.26. Графическое представление для решения краевой задачи

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

function res=MyBVP2(x,y) res=[y(2);y(2)];

end

Код достаточно простой и, думается, особых комментариев не требует. На рис. 6.27 этот же код представлен в окне редактора m-файлов.

Рис. 6.27. Функция для определения системы уравнений

267

Самоучитель Matlab

Функция для определения граничных условий такая:

function res=MyBC2(ya,yb) res=[ya(1)+1;yb(2)-yb(1)-2]; end

Код функции в окне редактора m-файлов представлен на рис. 6.28.

Рис. 6.28. Функция для определения граничных условий

В обоих случаях функции возвращают в качестве результата векторстолбец. Кроме этого, для формирования начального распределения для искомых функций создаем специальную функцию START(), которая в качестве результата возвращает вектор-строку с двумя элементами:

function res=START(x) res=[x.^2/2+x-1,x+1]; end

приближение для функции y(x)). Второй элемент соответствует функции y2(x) = x +1 (начальное приближение для производной y′(x) = x +1). Окно редактора m-файлов с кодом этой функции можно наблюдать на рис. 6.29.

Рис. 6.29. Функция для определения начального распределения

268

Глава 6. Интегрирование и дифференциальные уравнения

Рис. 6.30. Команды для решения краевой задачи

Рис. 6.31. Графическое представление для полученного решения

Далее используем следующую последовательность команд для решения краевой задачи:

>>IC2=bvpinit(0:0.01:1,@START);

>>sol=bvp5c(@MyBVP2,@MyBC2,IC2);

>>plot(sol.x,sol.y(1,:),'r-','LineWidth',2)

>>grid on

>>hold on

>>plot(sol.x,exp(sol.x)-2,'b:','LineWidth',3)

Инструкцией IC2=bvpinit(0:0.01:1,@START) создается начальное приближение для искомых функций. Первым аргументом функции bvpinit()

269

Самоучитель Matlab

передается список, формируемый командой 0:0.01:1. Вторым аргументом передается указатель @START на функцию, которой вычисляется начальное распределение для искомых функций. Далее это начальное условие, вместе с указателями @MyBVP2 и @MyBC2 на функцию уравнения и граничных условий, используется в команде sol=bvp5c(@MyBVP2,@MyBC2,IC2) для решения краевой задачи. Окно с использованными командами показано на рис. 6.30.

Как обычно, полученный результат отображается графически вместе с кривой для аналитического решения (рис. 6.31).

Традиционно кривая, построенная на основе числового решения, совпадает с кривой, построенной на основе известной зависимости для точного (аналитического) решения.

Завершающий пример

Да не было ничего! Все это происки!

В. Черномырдин

Еще один простой, но полезный пример использования встроенных функций для решения дифференциальных уравнений рассмотрим в этом разделе. Речь идет о том, что при определении некоторой функции используется функция для решения дифференциального уравнения. Другими словами, функция определяется так, что при вычислении результата решается дифференциальное уравнение. Обратимся к нижеприведенному коду:

>>f=@(x,y)(2*y./x+2*x.^3);

>>F=@(x,a,ya)(ode45(f,a:(x-a)/100:x,ya));

>>[x,y]=F(3,1,0);

>>plot(x,y,'r-','LineWidth',2)

>>grid on

>>hold on

>>plot(x,x.^2.*(x.^2-1),'b:','LineWidth',3)

Командой f=@(x,y)(2*y./x+2*x.^3) определяется функция, задающая

правую часть дифференциального уравнения y′(x) = 2 y(x) + 2x3 . Общее x

решение уравнения имеет вид y(x) = x2(x2 +C). Константа C вычисляется исходя из начальных условий.

Командой F=@(x,a,ya)(ode45(f,a:(x-a)/100:x,ya)) определяется функция, в которой вычисляется решение для упомянутого выше уравнения с начальным значением ya (третий аргумент функции F), которое вычисляется в точке a (второй аргумент функции F). Значения для решения дифференциального уравнения вычисляются от начального значения независимой переменной до того, что передано первым аргументом функции F.

270