1 <?xml version="1.0" encoding="UTF-8"?>
3 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
4 * Copyright (C) 2008 - INRIA
7 * This file must be used under the terms of the CeCILL.
8 * This source file is licensed as described in the file COPYING, which
9 * you should have received as part of this distribution. The terms
10 * are also available at
11 * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
14 <refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:svg="http://www.w3.org/2000/svg" xmlns:ns4="http://www.w3.org/1999/xhtml" xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org" xml:id="bvode" xml:lang="ru">
16 <refname>bvode</refname>
17 <refpurpose>задачи граничных значений для ОДУ с помощью метода коллокации
20 <refnamediv xml:id="bvodeS">
21 <refname>bvodeS</refname>
22 <refpurpose>упрощённый вызов bvode</refpurpose>
25 <title>Последовательность вызова</title>
27 zu=bvode(xpoints,N,m,x_low,x_up,zeta,ipar,ltol,tol,fixpnt,fsub,dfsub,gsub,dgsub,guess)
30 zu=bvodeS(xpoints,m,N,x_low,x_up,fsub,gsub,zeta, <optional_args>)
34 <title>Аргументы</title>
40 вектор-столбец длиной <literal>M</literal>. Решение ОДУ, вычисленное по сетке,
41 заданной точками. Она содержит <literal>z(u(x))</literal> для каждой искомой точки;
49 массив, который содержит точки для которых нужно найти решение;
57 скаляр с целочисленным значением, количество дифференциальных уравнений
58 (<literal>N</literal> <= 20).
66 вектор размером <literal>N</literal> с целочисленными элементами. Это вектор порядка
67 каждого из дифференциальных уравнений: <literal>m(i)</literal> указывает порядок
68 <literal>i</literal>-го дифференциального уравнения. Далее, <literal>M</literal>
69 будет представлять сумму элементов <literal>m</literal>. <!--Порядок
70 <literal>i</literal>-го дифференциального уравнения <literal>m(i)</literal> не
71 может быть больше 4.-->
78 <para>скаляр: левый конец интервала</para>
84 <para>скаляр: правый конец интервала</para>
91 вектор размером <literal>M</literal>; <literal>zeta(j)</literal> указывает
92 <literal>j</literal>-тую точку граничного условия (граничную точку). Необходимо,
93 чтобы выполнялось условие
94 <literal>x_low<=zeta(j)<=zeta(j+1)<=x_up</literal>.
97 Все точки граничных условий должны составлять сетку точек во всех используемых
98 сетках, см. ниже описание <literal>ipar(11)</literal> и <literal>fixpnt</literal>.
105 <para>массив с 11-ю целочисленными элементами:</para>
107 <literal>[nonlin, collpnt, subint, ntol, ndimf, ndimi, iprint,
108 iread, iguess, rstart,nfxpnt]
113 <term>nonlin: ipar(1)</term>
115 <para>0, если задача линейна; 1, если задача нелинейна.
120 <term>collpnt: ipar(2)</term>
123 Задаёт количество рядом расположенных точек на подынтервале, где
124 <literal>max(m(j)) <= collpnt <= 7</literal>.
127 Если <literal>ipar(2)=0</literal>, то <literal>collpnt</literal> установлен равным <literal>max( max(m(j))+1, 5-max(m(j)) )</literal>.
132 <term>subint: ipar(3)</term>
135 Задаёт количество подынтервалов в исходной сетке. Если
136 <literal>ipar(3) = 0</literal>, то <literal>bvode</literal> произвольным образом устанавливается <literal>subint = 5</literal>.
141 <term>ntol: ipar(4)</term>
144 Задаёт количество решений и допустимые отклонения производных.
145 Требуется, чтобы <literal>0 < ntol <= M</literal>.
146 <literal>ipar(4)</literal> должно быть установлено равным размеру аргумента
147 <literal>tol</literal> или равным <literal>0</literal>. В последнем случае
148 фактическое значение будет автоматически установлено равным
149 <literal>size(tol,'*')</literal>.
154 <term>ndimf: ipar(5)</term>
157 Задаёт размер <literal>fspace</literal> (вещественнозначный рабочий массив). Его значение
158 указывает ограничение по максимальному количеству подынтервалов
159 <literal>nmax</literal>.
162 Значение <literal>ipar(5)</literal> должно соответствовать ограничению
163 <literal>ipar(5)>=nmax*nsizef</literal>, где
164 <literal>nsizef=4 + 3*M + (5+collpnt*N)*(collpnt*N+M) + (2*M-
167 (<literal>nrec</literal> - количество граничных условий
173 <term>ndimi: ipar(6)</term>
176 Задаёт размер <literal>ispace</literal> (целочисленный рабочий массив).
177 Это значение указывает ограничение по <literal>nmax</literal>, максимальному количеству подынтервалов.
180 Значение <literal>ipar(6)</literal> должно соответствовать ограничению
181 <literal>ipar(6)>=nmax*nsizei</literal>, где <literal>nsizei= 3 + collpnt*N + M</literal>.
186 <term>iprint: ipar(7)</term>
188 <para>контроль на выходе, может принимать следующие значения:</para>
193 <para>для полной диагностической распечатки</para>
199 <para>для избранной распечатки</para>
205 <para>для отказа от распечатки</para>
212 <term>iread: ipar(8)</term>
219 заставляет <literal>bvode</literal> генерировать равномерную исходную
228 другие значения в Scilab'е пока не реализованы
236 если исходная сетка указана пользователем, то она будет определена в
237 <literal>fspace</literal> следующим образом: сетка будет занимать
238 <literal>fspace(1), ..., fspace(n+1)</literal>. Пользователю нужно предоставить
239 только внутренние точки сетки <literal>fspace(j) = x(j), j = 2, ..., n</literal>.
247 если исходная сетка представлена пользователем как в случае
248 <literal>ipar(8)=1</literal>, и, к тому же, никакого выбора адаптивной сетки не
257 <term>iguess: ipar(9)</term>
264 если никакого первоначального предположения для решения не
273 если первоначальное предположение предоставлено пользователем с
274 помощью аргумента <literal>guess</literal>.
282 если исходная сетка и коэффициенты приближённого решения предоставлены пользователем в <literal>fspace</literal> (прежняя и новая сетки одни и те же).
290 если прежняя сетка и коэффициенты приближённого решения предоставлены
291 пользователем в <literal>fspace</literal>, а новая сетка взята в два раза реже, то есть каждая вторая точка из прежней сетки.
299 если в дополнение к прежней исходной сетке и коэффициента приближённого
300 решения также предоставлена новая сетка в <literal>fspace</literal>
301 (см. описание выходных данных для более подробной информации по iguess = 2,
310 <term>ireg: ipar(10)</term>
316 <para>если задача является регулярной</para>
323 если первый коэффициент релаксации равен <literal>ireg</literal>, и
324 нелинейная итерация не зависит от прошлой сходимости (использовать только
325 для сверхчувствительных нелинейных задач)
334 если вы хотите возврата немедленно при (а) двух, следующих друг за другом,
335 несходимостях, либо (б) после получения ошибочной оценки в первый раз.
343 <term>nfxpnt: ipar(11)</term>
346 указывает количество фиксированных точек в сетке, отличных от
347 <literal>x_low</literal> и <literal>x_up</literal> (размерность <literal>fixpnt</literal>).
348 <literal>ipar(11)</literal> должна быть установлена равной размерности аргумента
349 <literal>fixpnt</literal> или равной <literal>0</literal>. В последнем случае фактическое значение будет автоматически установлено равным
350 <literal>size(fixpnt,'*')</literal>.
361 массив размерности <literal>ntol=ipar(4)</literal>.
362 <literal>ltol(j) = l</literal> определяет, что <literal>j</literal>-тый допуск в массиве
363 <literal> tol</literal> управляет ошибкой в <literal>l</literal>-том элементе
371 <mml:mo mml:stretchy="false">(</mml:mo><mml:mi>u</mml:mi><mml:mo mml:stretchy="false">)</mml:mo>
379 Также требуется, чтобы:
382 <literal>1 <= ltol(1) < ltol(2) < ... < ltol(ntol)
392 массив размерности <literal>ntol=ipar(4)</literal>.
395 <literal>tol(j)</literal> допуск ошибки в
396 <literal>ltol(j)</literal>-том элементе
404 <mml:mo mml:stretchy="false">(</mml:mo><mml:mi>u</mml:mi><mml:mo mml:stretchy="false">)</mml:mo>
411 . Таким образом код пытается удовлетворить
419 <mml:mfenced mml:close="∣" mml:open="∣">
422 <mml:mo mml:stretchy="false">(</mml:mo>
427 <mml:mo mml:stretchy="false">(</mml:mo><mml:mi>v</mml:mi><mml:mo mml:stretchy="false">)</mml:mo>
429 <mml:mo mml:stretchy="false">−</mml:mo><mml:mi>z</mml:mi>
432 <mml:mo mml:stretchy="false">(</mml:mo><mml:mi>u</mml:mi><mml:mo mml:stretchy="false">)</mml:mo>
435 <mml:mo mml:stretchy="false">)</mml:mo>
438 <mml:mi mml:fontstyle="italic">ltol</mml:mi>
440 <mml:mo mml:stretchy="false">(</mml:mo><mml:mi>j</mml:mi><mml:mo mml:stretchy="false">)</mml:mo>
445 <mml:mo mml:stretchy="false">≤</mml:mo><mml:mi mml:fontstyle="italic">tol</mml:mi>
450 <mml:mo mml:stretchy="false">(</mml:mo><mml:mi>j</mml:mi><mml:mo mml:stretchy="false">)</mml:mo>
452 <mml:mo mml:stretchy="false">⋅</mml:mo>
453 <mml:mfenced mml:close="∣" mml:open="∣">
456 <mml:mo mml:stretchy="false">(</mml:mo>
460 <mml:mo mml:stretchy="false">(</mml:mo><mml:mi>u</mml:mi><mml:mo mml:stretchy="false">)</mml:mo>
463 <mml:mo mml:stretchy="false">)</mml:mo>
466 <mml:mi mml:fontstyle="italic">ltol</mml:mi>
468 <mml:mo mml:stretchy="false">(</mml:mo><mml:mi>j</mml:mi><mml:mo mml:stretchy="false">)</mml:mo>
474 <mml:mo mml:stretchy="false">+</mml:mo><mml:mi mml:fontstyle="italic">tol</mml:mi>
477 <mml:mo mml:stretchy="false">(</mml:mo><mml:mi>j</mml:mi><mml:mo mml:stretchy="false">)</mml:mo>
479 <mml:mi>,</mml:mi><mml:mi mml:fontstyle="normal">for</mml:mi>
481 <mml:mi>j</mml:mi><mml:mo mml:stretchy="false">=</mml:mo><mml:mn>1</mml:mn>
483 <mml:mi mml:fontstyle="normal">:</mml:mi><mml:mi mml:fontstyle="normal">ntol</mml:mi>
485 <mml:annotation mml:encoding="StarMath 5.0">abs(z(v)-z(u))_{ltol(j)}
486 <= tol(j) cdot abs(z(u))_{ltol(j)} + tol(j), for
494 в каждом подынтервале, где
508 - вектор приближённого решения, а
520 - точное решение (неизвестное).
528 массив размером <literal>nfxpnt=ipar(11)</literal>. Он содержит точки отличные от
529 <literal>x_low</literal> и <literal>x_up</literal>, которые нужно включить во все
530 сетки. Код требует, чтобы все точки дополнительных условий, отличные от
531 <literal>x_low</literal> и <literal>x_up</literal> (см. описание <literal>zeta</literal> )
532 были включены в качестве фиксированных точек в <literal>fixpnt</literal>.
540 <link linkend="external">Внешняя функция</link>, используемая для вычисления вектор-столбца <literal>f=</literal>
547 <mml:mo mml:stretchy="false">[</mml:mo>
550 <mml:mi>f</mml:mi><mml:mn>1</mml:mn>
553 <mml:mo mml:stretchy="false">(</mml:mo>
555 <mml:mi>x</mml:mi><mml:mi>,</mml:mi><mml:mi>z</mml:mi>
557 <mml:mo mml:stretchy="false">(</mml:mo>
561 <mml:mo mml:stretchy="false">(</mml:mo><mml:mi>x</mml:mi><mml:mo mml:stretchy="false">)</mml:mo>
564 <mml:mo mml:stretchy="false">)</mml:mo>
567 <mml:mo mml:stretchy="false">)</mml:mo>
571 <mml:mi>f</mml:mi><mml:mn>2</mml:mn>
574 <mml:mo mml:stretchy="false">(</mml:mo>
576 <mml:mi>x</mml:mi><mml:mi>,</mml:mi><mml:mi>z</mml:mi>
578 <mml:mo mml:stretchy="false">(</mml:mo>
582 <mml:mo mml:stretchy="false">(</mml:mo><mml:mi>x</mml:mi><mml:mo mml:stretchy="false">)</mml:mo>
585 <mml:mo mml:stretchy="false">)</mml:mo>
588 <mml:mo mml:stretchy="false">)</mml:mo>
590 <mml:mi>,</mml:mi><mml:mo mml:stretchy="false">⋯</mml:mo><mml:mi>,</mml:mi>
592 <mml:mi>f</mml:mi><mml:mi>N</mml:mi>
595 <mml:mo mml:stretchy="false">(</mml:mo>
597 <mml:mi>x</mml:mi><mml:mi>,</mml:mi><mml:mi>z</mml:mi>
599 <mml:mo mml:stretchy="false">(</mml:mo>
603 <mml:mo mml:stretchy="false">(</mml:mo><mml:mi>x</mml:mi><mml:mo mml:stretchy="false">)</mml:mo>
606 <mml:mo mml:stretchy="false">)</mml:mo>
609 <mml:mo mml:stretchy="false">)</mml:mo>
612 <mml:mo mml:stretchy="false">]</mml:mo>
614 <mml:annotation mml:encoding="StarMath 5.0">[f_1(x,z(u(x)));
615 f_2(x,z(u(x))),dotsaxis,
624 для всех <literal>x</literal> таких, что <literal>x_low</literal> <= <literal>x</literal> <= <literal>x_up</literal>
625 и для любых <literal>z=z(u(x))</literal> (см. описание ниже).
627 <para>Внешняя функция должна иметь заголовки:</para>
630 <para>В Fortran последовательность вызова должна быть:</para>
631 <programlisting role="no-scilab-exec"><![CDATA[
632 subroutine fsub(x,zu,f)
633 double precision zu(*), f(*),x
637 <para>В C прототип функции должен быть:</para>
638 <programlisting role="no-scilab-exec"><![CDATA[
639 void fsub(double *x, double *zu, double *f)
643 <para>А в Scilab'е:</para>
644 <programlisting role="no-scilab-exec"><![CDATA[
645 function f=fsub(x,zu,parameters)
655 <link linkend="external">Внешняя функция</link>, используемая для вычисления
656 якобиана от <literal>f(x,z(u))</literal> в точке <literal>x</literal>, где
657 <literal>z(u(x))</literal> определена как для <literal>fsub</literal>, и массив
658 <literal>df</literal> размером <literal>N</literal> на <literal>M</literal> должен быть заполнен частными производными от <literal>f</literal>:
663 <imagedata align="center">
664 <mml:math scilab:localized="true">
667 <mml:mi mml:fontstyle="italic">df</mml:mi>
670 <mml:mo mml:stretchy="false">(</mml:mo>
676 <mml:mo mml:stretchy="false">)</mml:mo>
678 <mml:mo mml:stretchy="false">=</mml:mo>
681 <mml:mo mml:stretchy="false">∂</mml:mo>
688 <mml:mo mml:stretchy="false">∂</mml:mo>
699 <mml:mo mml:stretchy="true">{</mml:mo>
705 <mml:mo mml:stretchy="false">=</mml:mo>
708 <mml:mi mml:fontstyle="normal">:</mml:mi>
709 <mml:mi mml:fontstyle="italic">N</mml:mi>
716 <mml:mo mml:stretchy="false">=</mml:mo>
719 <mml:mi mml:fontstyle="normal">:</mml:mi>
720 <mml:mi mml:fontstyle="italic">M</mml:mi>
726 <mml:annotation mml:encoding="StarMath 5.0">df(i,j)= {partial{f_i}} over
727 {partial{z_j}} ~ for ~ left lbrace binom{i=1:N}{j=1:M} right
736 <para>Внешняя функция должна иметь заголовки:</para>
739 <para>В Fortran вызывающая последовательность должна быть:</para>
740 <programlisting role="no-scilab-exec"><![CDATA[
741 subroutine dfsub(x,zu,df)
742 double precision zu(*), df(*),x
746 <para>В C прототип функции должен быть:</para>
747 <programlisting role="no-scilab-exec"><![CDATA[
748 void dfsub(double *x, double *zu, double *df)
752 <para>И в Scilab'е:</para>
753 <programlisting role="no-scilab-exec"><![CDATA[
754 function df=dfsub(x,zu,parameters)
764 <link linkend="external">Внешняя функция</link>, используемая для вычисления
772 <mml:mi>g</mml:mi><mml:mi>i</mml:mi>
775 <mml:mo mml:stretchy="false">(</mml:mo>
778 <mml:mo mml:stretchy="false">ζ</mml:mo><mml:mi>i</mml:mi>
780 <mml:mi>,</mml:mi><mml:mi>z</mml:mi>
782 <mml:mo mml:stretchy="false">(</mml:mo>
786 <mml:mo mml:stretchy="false">(</mml:mo>
788 <mml:mo mml:stretchy="false">ζ</mml:mo><mml:mi>i</mml:mi>
790 <mml:mo mml:stretchy="false">)</mml:mo>
793 <mml:mo mml:stretchy="false">)</mml:mo>
796 <mml:mo mml:stretchy="false">)</mml:mo>
799 <mml:annotation mml:encoding="StarMath 5.0">g_i(%zeta_i,z(u(%zeta_i)))</mml:annotation>
814 <mml:mo mml:stretchy="false">(</mml:mo>
818 <mml:mo mml:stretchy="false">(</mml:mo>
820 <mml:mo mml:stretchy="false">ζ</mml:mo><mml:mi>i</mml:mi>
822 <mml:mo mml:stretchy="false">)</mml:mo>
825 <mml:mo mml:stretchy="false">)</mml:mo>
828 <mml:annotation mml:encoding="StarMath 5.0">z(u(%zeta_i))</mml:annotation>
834 <literal>z = zeta(i)</literal> для
835 <literal>1<=i<=M.</literal>
837 <para>Внешняя функция должна иметь заголовки:</para>
840 <para>В Fortran вызывающая последовательность должна быть:</para>
841 <programlisting role="no-scilab-exec"><![CDATA[
842 subroutine gsub(i,zu,g)
843 double precision zu(*), g(*)
848 <para>В C прототип функции должен быть:</para>
849 <programlisting role="no-scilab-exec"><![CDATA[
850 void gsub(int *i, double *zu, double *g)
854 <para>А в Scilab'е:</para>
855 <programlisting role="no-scilab-exec"><![CDATA[
856 function g=gsub(i,zu,parameters)
859 Заметьте, что в отличие от <literal>f</literal> в <literal>fsub</literal>, здесь только одно значение за вызов возвращается в <literal>g</literal>.
869 <link linkend="external">Внешняя функция</link>, используемая для вычисления
870 <literal>i</literal>-той строки якобиана от <literal>g(x,u(x))</literal>, где
871 <literal>z(u)</literal> такая, как для <literal>fsub</literal>, <literal>i</literal> как для <literal>gsub</literal> а <literal>M</literal>-вектор <literal>dg</literal> должен быть заполнен частными производными от <literal>g</literal>, то есть, для отдельного вызова вычисляется
876 <imagedata align="center">
877 <mml:math scilab:localized="true">
880 <mml:mi mml:fontstyle="italic">dg</mml:mi>
883 <mml:mo mml:stretchy="false">(</mml:mo>
889 <mml:mo mml:stretchy="false">)</mml:mo>
891 <mml:mo mml:stretchy="false">=</mml:mo>
894 <mml:mo mml:stretchy="false">∂</mml:mo>
901 <mml:mo mml:stretchy="false">∂</mml:mo>
912 <mml:mo mml:stretchy="true">{</mml:mo>
918 <mml:mo mml:stretchy="false">=</mml:mo>
921 <mml:mi mml:fontstyle="normal">:</mml:mi>
922 <mml:mi mml:fontstyle="italic">M</mml:mi>
929 <mml:mo mml:stretchy="false">=</mml:mo>
932 <mml:mi mml:fontstyle="normal">:</mml:mi>
933 <mml:mi mml:fontstyle="italic">M</mml:mi>
939 <mml:annotation mml:encoding="StarMath 5.0">dg(i,j)= {partial{g_i}} over
940 {partial{z_j}} ~ for ~ left lbrace binom{i=1:M}{j=1:M} right
949 <para>Внешняя функция должна иметь заголовки:</para>
952 <para>В Fortran вызывающая последовательность должна быть:</para>
953 <programlisting role="no-scilab-exec"><![CDATA[
954 subroutine dgsub(i,zu,dg)
955 double precision zu(*), dg(*)
959 <para>В C прототип функции должен быть:</para>
960 <programlisting role="no-scilab-exec"><![CDATA[
961 void dgsub(int *i, double *zu, double *dg)
965 <para>А в Scilab'е: </para>
966 <programlisting role="no-scilab-exec"><![CDATA[
967 function dg=dgsub(i,zu,parameters)
977 <link linkend="external">Внешняя функция</link>, используемая для вычисления
978 исходной аппроксимации для <literal>z(u(x))</literal> и
979 <literal>dmval(u(x))</literal>, вектора <literal>mj</literal>-тых производных от
980 <literal>u(x)</literal>. Заметьте, что эта процедура используется только если
981 <literal>ipar(9) = 1</literal>, и тогда все <literal>M</literal>
982 элементы <literal>zu</literal> и <literal>N</literal> элементы <literal>dmval</literal> должны быть вычислены для любого <literal>x</literal> такого, что <literal>x_low</literal> <= <literal>x</literal> <= <literal>x_up</literal>.
984 <para>Внешняя функция должна иметь заголовки:</para>
987 <para>В Fortran вызывающая последовательность должна быть:</para>
988 <programlisting role="no-scilab-exec"><![CDATA[
989 subroutine guess(x,zu,dmval)
990 double precision x,z(*), dmval(*)
994 <para>В C прототип функции должне быть:</para>
995 <programlisting role="no-scilab-exec"><![CDATA[
996 void fsub(double *x, double *zu, double *dmval)
1000 <para>А в Scilab'е: </para>
1001 <programlisting role="no-scilab-exec"><![CDATA[
1002 function [dmval,zu]=fsub(x,parameters)
1003 ]]></programlisting>
1009 <term><optional_args></term>
1011 <para>Необязательные аргументы, должны быть либо:</para>
1014 <para>любой левой частью упорядоченной последовательности значений:
1015 <literal>guess, dfsub, dgsub, fixpnt, ndimf, ndimi, ltol, tol,
1016 ntol,nonlin, collpnt, subint, iprint, ireg, ifail
1022 либо любой последовательностью <literal>arg_name=argvalue</literal>
1023 с <literal>arg_name</literal> из: <literal>guess</literal>,
1024 <literal>dfsub</literal>, <literal>dgsub</literal>,
1025 <literal>fixpnt</literal>, <literal>ndimf</literal>,
1026 <literal>ndimi</literal>, <literal>ltol</literal>,
1027 <literal>tol</literal>, <literal>ntol</literal>,
1028 <literal>nonlin</literal>, <literal>collpnt</literal>,
1029 <literal>subint</literal>, <literal>iprint</literal>,
1030 <literal>ireg</literal>, <literal>ifail</literal>.
1035 Все эти аргументы за исключением <literal>ifail</literal> описаны выше.
1036 <literal>ifail</literal> может использоваться для отображения вызова
1037 bvode в соответствии с выбранными необязательными аргументами. Если
1038 <literal>guess</literal>задано, то <literal>iguess</literal> равно 1.
1045 <title>Описание</title>
1047 Эти функции решают задачу многоточечных граничных значений для системы ОДУ смешанного порядка, заданной как:
1052 <imagedata align="left">
1056 <mml:mo mml:stretchy="true">{</mml:mo>
1065 <mml:mo mml:stretchy="false">(</mml:mo>
1070 <mml:mo mml:stretchy="false">)</mml:mo>
1073 <mml:mo mml:stretchy="false">=</mml:mo>
1080 <mml:mo mml:stretchy="false">(</mml:mo>
1086 <mml:mo mml:stretchy="false">(</mml:mo>
1090 <mml:mo mml:stretchy="false">(</mml:mo>
1092 <mml:mo mml:stretchy="false">)</mml:mo>
1095 <mml:mo mml:stretchy="false">)</mml:mo>
1098 <mml:mo mml:stretchy="false">)</mml:mo>
1102 <mml:mo mml:stretchy="true">{</mml:mo>
1108 <mml:mo mml:stretchy="false">=</mml:mo>
1111 <mml:mi mml:fontstyle="normal">:</mml:mi>
1118 <mml:mo mml:stretchy="false">∈</mml:mo>
1119 <mml:mfenced mml:close="[" mml:open="]">
1149 <mml:mo mml:stretchy="false">(</mml:mo>
1152 <mml:mo mml:stretchy="false">ζ</mml:mo>
1158 <mml:mo mml:stretchy="false">(</mml:mo>
1162 <mml:mo mml:stretchy="false">(</mml:mo>
1164 <mml:mo mml:stretchy="false">ζ</mml:mo>
1167 <mml:mo mml:stretchy="false">)</mml:mo>
1170 <mml:mo mml:stretchy="false">)</mml:mo>
1173 <mml:mo mml:stretchy="false">)</mml:mo>
1175 <mml:mo mml:stretchy="false">=</mml:mo>
1180 <mml:mtext>,</mml:mtext>
1182 <mml:mo mml:stretchy="false">=</mml:mo>
1185 <mml:mi mml:fontstyle="normal">:</mml:mi>
1191 <mml:annotation mml:encoding="StarMath 5.0">left lbrace
1192 stack{u_i^(m_i)=f_i(x,z(u(x))) ~~left lbrace { binom{i=1:N}{x
1193 in left ] a_l, a_u right [} } right none# `#
1194 g_j(%zeta_j,z(u(%zeta_j))) = 0~~j=1:M} right
1206 <imagedata align="left">
1210 <mml:mo mml:stretchy="true">{</mml:mo>
1215 <mml:mo mml:stretchy="false">=</mml:mo>
1218 <mml:mo mml:stretchy="false">∑</mml:mo>
1221 <mml:mo mml:stretchy="false">=</mml:mo>
1241 <mml:mo mml:stretchy="false">(</mml:mo>
1245 <mml:mo mml:stretchy="false">(</mml:mo>
1247 <mml:mo mml:stretchy="false">)</mml:mo>
1250 <mml:mo mml:stretchy="false">)</mml:mo>
1252 <mml:mo mml:stretchy="false">=</mml:mo>
1254 <mml:mo mml:stretchy="false">[</mml:mo>
1261 <mml:mo mml:stretchy="false">(</mml:mo>
1263 <mml:mo mml:stretchy="false">)</mml:mo>
1272 <mml:mo mml:stretchy="false">(</mml:mo>
1274 <mml:mo mml:stretchy="false">)</mml:mo>
1277 <mml:mo mml:stretchy="false">⋯</mml:mo>
1282 <mml:mo mml:stretchy="false">(</mml:mo>
1288 <mml:mo mml:stretchy="false">−</mml:mo>
1291 <mml:mo mml:stretchy="false">)</mml:mo>
1295 <mml:mo mml:stretchy="false">(</mml:mo>
1297 <mml:mo mml:stretchy="false">)</mml:mo>
1301 <mml:mo mml:stretchy="false">⋯</mml:mo>
1309 <mml:mo mml:stretchy="false">(</mml:mo>
1311 <mml:mo mml:stretchy="false">)</mml:mo>
1320 <mml:mo mml:stretchy="false">(</mml:mo>
1322 <mml:mo mml:stretchy="false">)</mml:mo>
1325 <mml:mo mml:stretchy="false">⋯</mml:mo>
1331 <mml:mo mml:stretchy="false">(</mml:mo>
1337 <mml:mo mml:stretchy="false">−</mml:mo>
1340 <mml:mo mml:stretchy="false">)</mml:mo>
1344 <mml:mo mml:stretchy="false">(</mml:mo>
1346 <mml:mo mml:stretchy="false">)</mml:mo>
1349 <mml:mo mml:stretchy="false">]</mml:mo>
1366 <mml:mo mml:stretchy="false">≤</mml:mo>
1368 <mml:mo mml:stretchy="false">ζ</mml:mo>
1372 <mml:mo mml:stretchy="false">≤</mml:mo>
1373 <mml:mo mml:stretchy="false">⋯</mml:mo>
1375 <mml:mo mml:stretchy="false">≤</mml:mo>
1377 <mml:mo mml:stretchy="false">ζ</mml:mo>
1381 <mml:mo mml:stretchy="false">≤</mml:mo>
1390 <mml:annotation mml:encoding="StarMath 5.0">left lbrace stack
1391 { M= sum from {i=1} to {N} m_i # ~# z(u(x)) =
1392 [u_1(x);u_1^{1}(x);dotsaxis u_1^(m_1-1)(x);~ dotsaxis;
1393 ~u_N(x);u_N^{1}(x);dotsaxis; u_N^(m_N-1)(x)]# ~# x_l <=
1394 %zeta_1 <= dotsaxis <= %zeta_M <=x_u}right
1403 Аргумент <literal>zu</literal>, используемый внешними функциями и возвращаемый
1404 <literal>bvode</literal>, является вектор-столбцом, сформированным элементами
1405 <literal>z(u(x))</literal> для заданных <literal>x</literal>.
1408 Метод, используемый для аппроксимации решения, <literal>u</literal> является
1409 коллокацией в гауссовских точках, требующих <literal>m(i)-1</literal> непрерывных
1410 производных в <literal>i</literal>-том элементе, <literal>i = 1:N</literal>. Здесь,
1411 <literal>k</literal> - количество точек коллокации (этапов) на подынтервале, и оно
1412 выбирается так, чтобы <literal>k .ge. max m(i)</literal>. Используется представление
1413 решения одночленного решения Рунге-Кутты.
1417 <title>Примеры</title>
1418 <para>Первые две задачи взяты из [1] раздела Литература.
1423 <emphasis role="bold">Задача 1</emphasis> описывает однородно нагруженную балку переменной жёсткости просто поддерживаемую на обоих концах.
1425 <para>Она может быть определена следующим образом.</para>
1426 <para>Решим дифференциальное уравнение четвёртого порядка:</para>
1429 <imagedata align="left">
1449 <mml:mo mml:stretchy="false">(</mml:mo>
1451 <mml:mo mml:stretchy="false">)</mml:mo>
1453 <mml:mo mml:stretchy="false">=</mml:mo>
1458 <mml:mo mml:stretchy="false">−</mml:mo>
1462 <mml:mo mml:stretchy="false">.</mml:mo>
1468 <mml:mo mml:stretchy="false">.</mml:mo>
1487 <mml:mo mml:stretchy="false">(</mml:mo>
1489 <mml:mo mml:stretchy="false">)</mml:mo>
1491 <mml:mo mml:stretchy="false">−</mml:mo>
1495 <mml:mo mml:stretchy="false">.</mml:mo>
1498 <mml:mo mml:stretchy="false">.</mml:mo>
1516 <mml:mo mml:stretchy="false">(</mml:mo>
1518 <mml:mo mml:stretchy="false">)</mml:mo>
1528 <mml:annotation mml:encoding="StarMath 5.0">d^4 over {d
1529 x^4} u(x)={1 -6*x^2*{{d^3} over {d x^3}} u(x)-6*x*{{d^2}
1530 over {d x^2}} u(x)} over x^3
1537 <para>Налагаемые граничные условия:</para>
1540 <imagedata align="left">
1544 <mml:mo mml:stretchy="true">{</mml:mo>
1551 <mml:mo mml:stretchy="false">(</mml:mo>
1553 <mml:mo mml:stretchy="false">)</mml:mo>
1555 <mml:mo mml:stretchy="false">=</mml:mo>
1577 <mml:mo mml:stretchy="false">(</mml:mo>
1579 <mml:mo mml:stretchy="false">)</mml:mo>
1583 <mml:mo mml:stretchy="false">(</mml:mo>
1585 <mml:mo mml:stretchy="false">)</mml:mo>
1587 <mml:mo mml:stretchy="false">=</mml:mo>
1597 <mml:mo mml:stretchy="false">(</mml:mo>
1599 <mml:mo mml:stretchy="false">)</mml:mo>
1601 <mml:mo mml:stretchy="false">=</mml:mo>
1623 <mml:mo mml:stretchy="false">(</mml:mo>
1625 <mml:mo mml:stretchy="false">)</mml:mo>
1629 <mml:mo mml:stretchy="false">(</mml:mo>
1631 <mml:mo mml:stretchy="false">)</mml:mo>
1633 <mml:mo mml:stretchy="false">=</mml:mo>
1640 <mml:annotation mml:encoding="StarMath 5.0">left lbrace
1641 stack{ u(1)=0# {{d^2} over {d x^2}} u(x)(1)=0# u(2)=0#
1642 {{d^2} over {d x^2}} u(x)(2)=0} right
1650 <para>Известно точное решение этой задачи:</para>
1653 <imagedata align="left">
1660 <mml:mo mml:stretchy="false">(</mml:mo>
1662 <mml:mo mml:stretchy="false">)</mml:mo>
1664 <mml:mo mml:stretchy="false">=</mml:mo>
1671 <mml:mo mml:stretchy="false">(</mml:mo>
1675 <mml:mi>log</mml:mi>
1678 <mml:mo mml:stretchy="false">(</mml:mo>
1680 <mml:mo mml:stretchy="false">)</mml:mo>
1682 <mml:mo mml:stretchy="false">−</mml:mo>
1686 <mml:mo mml:stretchy="false">)</mml:mo>
1691 <mml:mo mml:stretchy="false">(</mml:mo>
1694 <mml:mo mml:stretchy="false">−</mml:mo>
1697 <mml:mo mml:stretchy="false">)</mml:mo>
1699 <mml:mo mml:stretchy="false">+</mml:mo>
1706 <mml:mo mml:stretchy="false">[</mml:mo>
1712 <mml:mo mml:stretchy="false">+</mml:mo>
1715 <mml:mo mml:stretchy="false">(</mml:mo>
1718 <mml:mo mml:stretchy="false">+</mml:mo>
1721 <mml:mo mml:stretchy="false">)</mml:mo>
1724 <mml:mi>log</mml:mi>
1727 <mml:mo mml:stretchy="false">(</mml:mo>
1729 <mml:mo mml:stretchy="false">)</mml:mo>
1731 <mml:mo mml:stretchy="false">−</mml:mo>
1736 <mml:mo mml:stretchy="false">]</mml:mo>
1739 <mml:annotation mml:encoding="StarMath 5.0">u(x)=1 over 4
1740 (10`log(2)-3)`(1-x)+ {1 over 2} [ {1 over x} + {(3+x) `
1748 <programlisting role="example"><![CDATA[
1749 N=1;// только одно дифференциальное уравнение
1750 m=4;//дифференциальное уравнение четвёртого порядка
1754 x_up=2; // пределы x
1755 zeta=[x_low,x_low,x_up,x_up]; // два ограничения (на значение u и его вторую производную) на каждой границе.
1758 //Эти функции вызываются программой решения с zu=[u(x);u'(x);u''(x);u'''(x)]
1760 // - Функция, которая вычисляет правую сторону дифференциального уравнения
1761 function f=fsub(x,zu)
1762 f=(1-6*x^2*zu(4)-6*x*zu(3))/x^3
1765 // - Функция, которая вычисляет производную fsub по zu
1766 function df=dfsub(x,zu)
1767 df=[0,0,-6/x^2,-6/x]
1770 // - Функция, которая вычисляет i-тое ограничение для заданной i
1771 function g=gsub(i,zu),
1773 case 1 then //x=zeta(1)=1
1775 case 2 then //x=zeta(2)=1
1777 case 3 then //x=zeta(3)=2
1779 case 4 then //x=zeta(4)=2
1784 // - Функция, которая вычисляет производную gsub по z
1785 function dg=dgsub(i,z)
1787 case 1 then //x=zeta(1)=1
1789 case 2 then //x=zeta(2)=1
1791 case 3 then //x=zeta(3)=2
1793 case 4 then //x=zeta(4)=2
1798 // - Функция, которая вычисляет начальное предположение, здесь не используется
1799 function [zu,mpar]=guess(x)
1804 //определим функцию, которая вычисляет точное значение u для заданного x (в целях проверки)
1805 function zu=trusol(x)
1807 zu(1) = 0.25*(10*log(2)-3)*(1-x) + 0.5 *( 1/x + (3+x)*log(x) - x)
1808 zu(2) = -0.25*(10*log(2)-3) + 0.5 *(-1/x^2 + (3+x)/x + log(x) - 1)
1809 zu(3) = 0.5*( 2/x^3 + 1/x - 3/x^2)
1810 zu(4) = 0.5*(-6/x^4 - 1/x/x + 6/x^3)
1813 fixpnt=[ ];//Все граничные условия находятся в x_low и x_up
1815 // nonlin collpnt n ntol ndimf ndimi iprint iread iguess rstart nfxpnt
1816 +ipar=[1 collpnt 10 4 nmax*nsizef nmax*nsizei -1 0 0 0 0 ]
1818 ltol=[1,3];//установка контроля допуска на zu(1) и zu(3)
1819 tol=[1.e-11,1.e-11];//установка значений допуска для контроля допуска
1820 xpoints=x_low:0.01:x_up;
1822 zu=bvode(xpoints,N,m,x_low,x_up,zeta,ipar,ltol,tol,fixpnt,...
1823 fsub,dfsub,gsub,dgsub,guess)
1824 //проверка ограничений
1825 zu([1,3],[1 $]) // должно быть нулевым
1826 plot(xpoints,zu(1,:)) // процесс решения u
1829 zu1=[zu1,trusol(x)];
1832 ]]></programlisting>
1836 Та же задача с использованием <literal>bvodeS</literal> и начальным предположением.
1838 <programlisting role="no-scilab-exec"><![CDATA[
1839 function [z,lhS]=zstart(x)
1840 z=zeros(5,1);z(5)=1;
1843 zu=bvode(xpoints,N,m,x_low,x_up,zeta,ltol=ltol,tol=tol,guess=zstart)
1844 ]]></programlisting>
1848 <emphasis role="bold">Задача 2</emphasis> описывает малую конечную деформацию тонкой мелкой сферической чаши постоянной толщины, подверженной квадратично меняющемуся асимметричному распределению внешнего давления. Здесь <latex>$\varphi$</latex> описывает изменение угла по меридиану деформированной стенки, а <latex>$\psi$</latex> является функцией напряжений. Для <latex>$\varepsilon=\mu=10^{-3}$</latex> можно найти два различных решения в зависимости от точки начала.
1852 <imagedata align="left">
1856 <mml:mo mml:stretchy="true">{</mml:mo>
1862 <mml:mo mml:stretchy="false">ε</mml:mo>
1865 <mml:mo mml:stretchy="false">μ</mml:mo>
1869 <mml:mo mml:stretchy="false">(</mml:mo>
1871 <mml:mo mml:stretchy="false">φ</mml:mo>
1876 <mml:mo mml:stretchy="false">+</mml:mo>
1879 <mml:mo mml:stretchy="false">φ</mml:mo>
1885 <mml:mo mml:stretchy="false">−</mml:mo>
1887 <mml:mo mml:stretchy="false">φ</mml:mo>
1895 <mml:mo mml:stretchy="false">)</mml:mo>
1897 <mml:mo mml:stretchy="false">+</mml:mo>
1898 <mml:mo mml:stretchy="false">ψ</mml:mo>
1903 <mml:mo mml:stretchy="false">(</mml:mo>
1906 <mml:mo mml:stretchy="false">−</mml:mo>
1908 <mml:mo mml:stretchy="false">φ</mml:mo>
1912 <mml:mo mml:stretchy="false">)</mml:mo>
1914 <mml:mo mml:stretchy="false">−</mml:mo>
1915 <mml:mo mml:stretchy="false">φ</mml:mo>
1917 <mml:mo mml:stretchy="false">+</mml:mo>
1918 <mml:mo mml:stretchy="false">γ</mml:mo>
1923 <mml:mo mml:stretchy="false">(</mml:mo>
1926 <mml:mo mml:stretchy="false">−</mml:mo>
1935 <mml:mo mml:stretchy="false">)</mml:mo>
1937 <mml:mo mml:stretchy="false">=</mml:mo>
1944 <mml:mo mml:stretchy="false">μ</mml:mo>
1947 <mml:mo mml:stretchy="false">(</mml:mo>
1949 <mml:mo mml:stretchy="false">ψ</mml:mo>
1954 <mml:mo mml:stretchy="false">+</mml:mo>
1957 <mml:mo mml:stretchy="false">ψ</mml:mo>
1963 <mml:mo mml:stretchy="false">−</mml:mo>
1965 <mml:mo mml:stretchy="false">ψ</mml:mo>
1973 <mml:mo mml:stretchy="false">)</mml:mo>
1975 <mml:mo mml:stretchy="false">−</mml:mo>
1976 <mml:mo mml:stretchy="false">φ</mml:mo>
1980 <mml:mo mml:stretchy="false">(</mml:mo>
1983 <mml:mo mml:stretchy="false">−</mml:mo>
1985 <mml:mo mml:stretchy="false">φ</mml:mo>
1989 <mml:mo mml:stretchy="false">)</mml:mo>
1991 <mml:mo mml:stretchy="false">=</mml:mo>
2000 <mml:mo mml:stretchy="false"><</mml:mo>
2003 <mml:mo mml:stretchy="false"><</mml:mo>
2009 <mml:annotation mml:encoding="StarMath 5.0">left lbrace
2010 stack{{%epsilon^4} over %mu (%phi''+ {%phi'} over x
2011 -{%phi} over {x^2} ) +%psi (1- %phi over x ) -%phi+%gamma
2012 x(1-{x^2} over 2)=0# %mu(%psi''+{%psi'} over x-%psi over
2013 {x^2})-%phi(1-%phi over{2x})=0 # 0<x<1} right
2021 <para>Граничные условия:</para>
2024 <imagedata align="left">
2025 <mml:math scilab:localized="true">
2028 <mml:mo mml:stretchy="true">{</mml:mo>
2032 <mml:mo mml:stretchy="false">φ</mml:mo>
2033 <mml:mo mml:stretchy="false">=</mml:mo>
2041 <mml:mo mml:stretchy="false">ψ</mml:mo>
2044 <mml:mo mml:stretchy="false">−</mml:mo>
2045 <mml:mn>0,3</mml:mn>
2049 <mml:mo mml:stretchy="false">ψ</mml:mo>
2050 <mml:mo mml:stretchy="false">+</mml:mo>
2051 <mml:mn>0,7x</mml:mn>
2053 <mml:mo mml:stretchy="false">=</mml:mo>
2060 <mml:annotation mml:encoding="StarMath 5.0">left lbrace
2061 stack{ %phi=0# x`%psi'-0.3%psi+0.7x=0} right
2070 для <literal>x=0</literal> и <literal>x=1</literal>
2072 <programlisting role="example"><![CDATA[
2073 N=2;// два дифференциальных уравнения
2074 m=[2 2];// каждое дифференциальное уравнение второго порядка
2077 x_low=0;x_up=1; // пределы по x
2078 zeta=[x_low,x_low, x_up x_up]; //два ограничения на каждой границе
2081 //Эти функции вызываются программой решения с zu=[u1(x);u1'(x);u2(x);u2'(x)]
2083 // - Функция, которая вычисляет правую сторону дифференциального уравнения
2084 function f=fsub2(x,zu,eps,dmu,eps4mu,gam,xt),
2085 f=[zu(1)/x^2-zu(2)/x+(zu(1)-zu(3)*(1-zu(1)/x)-gam*x*(1-x^2/2))/eps4mu //phi''
2086 zu(3)/x^2-zu(4)/x+zu(1)*(1-zu(1)/(2*x))/dmu];//psi''
2089 // - Функция, которая вычисляет производную от fsub по zu
2090 function df=dfsub2(x,zu,eps,dmu,eps4mu,gam,xt),
2091 df=[1/x^2+(1+zu(3)/x)/eps4mu, -1/x, -(1-zu(1)/x)/eps4mu, 0
2092 (1-zu(1)/x)/dmu 0 1/x^2 -1/x];
2095 // - Функция, которая вычисляет i-тое ограничение для заданного i
2096 function g=gsub2(i,zu),
2098 case 1 then //x=zeta(1)=0
2100 case 2 then //x=zeta(2)=0
2101 g=-0.3*zu(3) //x*psi'-0.3*psi+0.7x=0
2102 case 3 then //x=zeta(3)=1
2104 case 4 then //x=zeta(4)=1
2105 g=1*zu(4)-0.3*zu(3)+0.7*1 //x*psi'-0.3*psi+0.7x=0
2109 // - Функция, которая вычисляет производную от gsub по z
2110 function dg=dgsub2(i,z)
2112 case 1 then //x=zeta(1)=1
2114 case 2 then //x=zeta(2)=1
2116 case 3 then //x=zeta(3)=2
2118 case 4 then //x=zeta(4)=2
2127 xt=sqrt(2*(gam-1)/gam)
2129 fixpnt=[ ];// все граничные условия размещены в x_low и x_up
2131 nsizef=4+3*M+(5+collpnt*N)*(collpnt*N+M)+(2*M-2)*2*M ;
2132 nsizei=3 + collpnt*N+M;;
2134 // nonlin collpnt n ntol ndimf ndimi iprint iread iguess rstart nfxpnt
2135 ipar=[1 k 10 4 nmax*nsizef nmax*nsizei -1 0 0 0 0 ]
2137 ltol=1:4;//установка контроля допусков на zu(1), zu(2), zu(3) и zu(4)
2138 tol=[1.e-5,1.e-5,1.e-5,1.e-5];//установка значений для контроля допусков
2139 xpoints=x_low:0.01:x_up;
2141 // - Функция, которая вычисляет начальное предположение, здесь не используется
2142 function [zu,dmval]=guess2(x,gam),
2143 cons=gam*x*(1-x^2/2)
2144 dcons=gam*(1-3*x^2/2)
2148 zu=[0 0 -cons -dcons]
2151 zu=[2*x;2;-2*x+cons;-2*dcons]
2156 zu=bvode(xpoints,N,m,x_low,x_up,zeta,ipar,ltol,tol,fixpnt,...
2157 fsub2,dfsub2,gsub2,dgsub2,guess2);
2158 scf(1);clf();plot(xpoints,zu([1 3],:)) // процесс решения phi и psi
2160 //использование начального предположения
2163 zu2=bvode(xpoints,N,m,x_low,x_up,zeta,ipar,ltol,tol,fixpnt,...
2164 fsub2,dfsub2,gsub2,dgsub2,guess2);
2165 scf(2);clf();plot(xpoints,zu2([1 3],:)) // процесс решения phi и psi
2166 ]]></programlisting>
2170 <emphasis role="bold">Задача нахождения собственных значений</emphasis>
2172 <programlisting role="example"><![CDATA[
2174 // BV: y(0)=y'(0); y(1)=0
2175 // Собственные функции и собственные значения y(x,n)=sin(s(n)*(1-x)), la(n)=s(n)^2,
2176 // где s(n) - нули f(s,n)=s+atan(s)-(n+1)*pi, n=0,1,2,...
2177 // Чтобы получить третье граничное условие, мы выбираем y(0)=1
2178 // (с y(x) также c*y(x) является решением для каждой константы c.)
2179 // Решим следующую систему ОДУ:
2182 // BV: y(0)=y'(0), y(0)=1; y(1)=0
2183 // z=[y(x) ; y'(x) ; la]
2185 function rhs=fsub(x,z)
2189 function g=gsub(i,z)
2190 g=[z(1)-z(2) z(1)-1 z(1)]
2194 // Следующая начальная функция годится для первых 8 собственных функций.
2195 function [z,lhs]=ystart(x,z,la0)
2207 //Имеем s(n)-(n+1/2)*pi -> 0 для n, стремящимся к бесконечности.
2208 la0=input('n-ное собственное значение: n= ?');la0=(%pi/2+la0*%pi)^2;
2210 z=bvodeS(x,m,n,a,b,fsub,gsub,zeta,ystart=list(ystart,la0));
2211 // Тот же вызов без вывода на экран
2212 z=bvodeS(x,m,n,a,b,fsub,gsub,zeta,ystart=list(ystart,la0),iprint=1);
2213 // Тот же вызов с подробным выводом на экран
2214 z=bvodeS(x,m,n,a,b,fsub,gsub,zeta,ystart=list(ystart,la0),iprint=-1);
2217 plot(x,[z(1,:)' z(2,:)'])
2218 xtitle(['Начальное значение = '+string(la0);'Собственное значение = '+string(z(3,1))],'x',' ')
2219 legend(['y(x)';'y''(x)'])
2220 ]]></programlisting>
2221 <scilab:image localized="true">
2222 function rhs=fsub(x,z)
2226 function g=gsub(i,z)
2227 g=[z(1)-z(2) z(1)-1 z(1)]
2231 function [z,lhs]=ystart(x,z,la0)
2244 la0=(%pi/2+la0*%pi)^2;
2246 z=bvodeS(x,m,n,a,b,fsub,gsub,zeta,ystart=list(ystart,la0),iprint=1);
2248 plot(x,[z(1,:)' z(2,:)']);
2249 xtitle(['Начальное значение = '+string(la0);'Собственное значение = '+string(z(3,1))],'x',' ');
2250 legend(['y(x)';'y''(x)']);
2255 <emphasis role="bold">Краевая задача, у которой решений больше одного</emphasis>
2257 <programlisting role="example"><![CDATA[
2258 // DE: y''(x)=-exp(y(x))
2259 // BV: y(0)=0; y(1)=0
2260 // Эта краевая задача имеет более одного решения.
2261 // Показано как найти два из них с помощью некоторой
2262 // предварительной информации о решениях y(x) для построения функции ystart.
2274 function rhs=fsub(x,z),rhs=-exp(z(1));endfunction
2276 function g=gsub(i,z)
2281 function [z,lhs]=ystart(x,z,M)
2282 //z=[4*x*(1-x)*M ; 4*(1-2*x)*M]
2284 //lhs=[-exp(4*x*(1-x)*M)]
2290 z=bvodeS(x,m,n,a,b,fsub,gsub,zeta,ystart=list(ystart,M),tol=tol);
2292 z1=bvodeS(x,m,n,a,b,fsub,gsub,zeta,ystart=list(ystart,M),tol=tol);
2296 // Интегрирование ОДУ даёт, например, два решения yex и yex1.
2299 y=c.*(1-tanh(sqrt(c)/4).^2)-2;
2304 y=log(c/2*(1-tanh(sqrt(c)*(1/4-x/2)).^2))
2308 y=2*c1^2+tanh(1/4/c1)^2-1;
2312 function y=yex1(x,c1)
2313 y=log((1-tanh((2*x-1)/4/c1).^2)/2/c1/c1)
2316 disp(norm(z(1,:)-yex(x)),'norm(yex(x)-z(1,:))= ')
2317 disp(norm(z1(1,:)-yex1(x)),'norm(yex1(x)-z1(1,:))= ')
2320 plot2d(x,z(1,:),style=[5])
2321 xtitle('Два различных решения','x',' ')
2323 plot2d(x,z1(1,:),style=[5])
2325 ]]></programlisting>
2326 <scilab:image localized="true">
2336 function rhs=fsub(x,z)
2340 function g=gsub(i,z)
2345 function [z,lhs]=ystart(x,z,M)
2346 //z=[4*x*(1-x)*M ; 4*(1-2*x)*M]
2348 //lhs=[-exp(4*x*(1-x)*M)]
2354 z=bvodeS(x,m,n,a,b,fsub,gsub,zeta,ystart=list(ystart,M),tol=tol);
2356 z1=bvodeS(x,m,n,a,b,fsub,gsub,zeta,ystart=list(ystart,M),tol=tol);
2360 // Интегрирование ОДУ даёт, например, два решения yex и yex1.
2363 y=c.*(1-tanh(sqrt(c)/4).^2)-2;
2368 y=log(c/2*(1-tanh(sqrt(c)*(1/4-x/2)).^2))
2372 y=2*c1^2+tanh(1/4/c1)^2-1
2376 function y=yex1(x,c1)
2377 y=log((1-tanh((2*x-1)/4/c1).^2)/2/c1/c1)
2380 disp(norm(z(1,:)-yex(x)),'norm(yex(x)-z(1,:))= ')
2381 disp(norm(z1(1,:)-yex1(x)),'norm(yex1(x)-z1(1,:))= ')
2384 plot2d(x,z(1,:),style=[5])
2385 xtitle('Два различных решения','x',' ')
2387 plot2d(x,z1(1,:),style=[5])
2393 <emphasis role="bold">Многоточечная краевая задача</emphasis>
2395 <programlisting role="example"><![CDATA[
2397 // z=[y(x);y'(x);y''(x)]
2398 // BV: y(-1)=2 y(1)=2
2399 // Дополнительное условие: y(0)=1
2402 // Точка дополнительного условия c должна быть включена в массив fixpnt.
2406 function rhs=fsub(x,z)
2410 function g=gsub(i,z)
2411 g=[z(1)-2 z(1)-1 z(1)-2]
2419 z=bvodeS(x,m,n,a,b,fsub,gsub,zeta,fixpnt=c);
2422 y=x.^3/6+x.^2-x./6+1
2425 disp(norm(yex(x)-z(1,:)),'norm(yex(x)-z(1,:))= ')
2426 ]]></programlisting>
2430 <emphasis role="bold">
2431 Квантовое уравнение Неймана с 2 "собственными значениями"
2432 (c_1 и c2). Используется продолжение.
2435 <programlisting role="example"><![CDATA[
2437 // Квантовое уравнение Неймана с 2 "собственными значениями" c_1 и c_2
2438 // (c_1=v-c_2-c_3, где v - параметр, используемый в продолжении)
2440 // diff(f,x,2) + (1/2)*(1/x + 1/(x-1) + 1/(x-y))*diff(f,x)
2441 // - (c_1/x + c_2/(x-1) + c_3/(x-y))* f(x) = 0
2442 // diff(c_2,x)=0, diff(c_3,x) = 0
2444 // и 4 "граничных" условия: diff(f,x)(a_k)=2*c_k*f(a_k) для
2445 // k=1,2,3, a_k=(0, 1 , y) и нормировка f(1) = 1
2447 // Вектор z - это z_1=f, z_2=diff(f,x), z_3=c_2 и z_4=c_3
2448 // Предположение выбирается так, чтобы иметь один узел в интервале
2449 // [0,1], f(x)=2*x-1 такая, что f(1)=1, c_2 и c_3 выбираются так,
2450 // чтобы исключить полюсы в дифференциальном уравнении в 1.0 и
2451 // y, z_3=1, z_4=1/(2*y-1)
2452 // См.: http://arxiv.org/pdf/hep-th/0407005
2458 eigens=zeros(3,40); // Для хранения результатов
2460 // Общая настройка bvode
2462 // Количество дифференциальных уравнений
2465 // Порядки уравнений
2468 // Нелинейная задача
2471 // Количество точек коллокации
2474 // Исходная сетка из 4 подынтервалов
2478 // Размер fspace, ispace, см. colnew.f для изменения размера
2482 // Средний вывод данных
2485 // Есть исходное приближение
2488 // fixpnt - это массив, содержащий все фиксированные точки в сетке, в
2489 // частности, "граничные" точки, за исключением aleft и aright, ipar[11]
2490 // её размер, здесь только одна внутренняя "граничная точка"
2494 // Допустимые отклонения для всех элементов z_1, z_2, z_3, z_4
2497 // Проверка допустимого отклонения для f и diff(f,x) и для c_2 и c_3
2498 ltol = [1, 2, 3, 4];
2499 tol = [1d-5, 1d-5, 1d-5, 1d-5];
2501 // Определение дифференциальных уравнений
2503 function [f]=fsub(x,z)
2504 f = [ -.5*(1/x+1/(x-1)+1/(x-y))*z(2) +...
2505 z(1) * ((v-z(3)-z(4))/x + z(3)/(x-1) + z(4)/(x-y)),...
2508 function [df] = dfsub(x,z)
2509 df = [(v-z(3)-z(4))/x + z(3)/(x-1) + z(4)/(x-y),...
2510 -.5*(1/x+1/(x-1)+1/(x-y)),z(1)/(x*(x-1)),z(1)*y/(x*(x-y));...
2514 // Граничные условия
2516 function [g]=gsub(i,z)
2518 case 1, g = z(2) - 2*z(1)*(v-z(3)-z(4))
2519 case 2, g = z(2) - 2*z(1)*z(3)
2521 case 4, g = z(2) - 2*z(1)*z(4)
2524 function [dg]=dgsub(i,z)
2526 case 1, dg = [-2*(v-z(3)-z(4)),1.,2*z(1),2*z(1)]
2527 case 2, dg = [-2*z(3),1.,-2*z(1),0]
2528 case 3, dg = [1,0,0,0]
2529 case 4, dg = [-2*z(4),1.,0,-2*z(1)]
2533 // Начало вычисления
2535 // Положения краевых условий, отсортированные
2536 zeta = [0.0d0, 1.0d0, 1.0d0, y];
2541 // Массив из 40 значений v, исследуемых продолжением, и массив из 202
2542 // точек, в которых следует вычислить функцию f.
2543 valv = [linspace(0,.9,10) logspace(0,2,30)];
2544 res = [linspace(0,.99,100) linspace(1,y,101)];
2546 // Собственные состояния характеризуются количеством узлов на интервале
2547 // [0,1] и на [1,y], здесь предполжение выбирает один узел (ноль) на
2548 // интервале [0,1] с линейной функцией f(x)=2*x-1 и константами c_2, c_3,
2549 // так что dmval=0. Заметьте, что вектор z имеет mstar = 4 элементов,
2550 // в то время как dmval имеет ncomp = 3 элемента.
2552 function [z,dmval]=guess(x)
2553 z=[2*x-1, 2., 1., 1/(2*y-1)]
2557 // Первое выполнение имеет ipar(9)=1 и использует предположение
2558 // Последующие выполнения имеют ipar(9)=3 и используют продолжение.
2559 // Это выполняется в плотно закрытом цикле, чтобы не нарушить стек
2563 sol=bvode(res,ncomp,m,aleft,aright,zeta,ipar,ltol,tol,fixpnt,...
2564 fsub,dfsub,gsub,dgsub,guess);
2565 eigens(:,i)=[v;sol(3,101);sol(4,101)]; // c_2 и c_3 являются константами!
2570 // Чтобы увидеть изменение собственных значений с v, disp(eigens)
2571 // Обратите внимание на то, что они изменяются гладко.
2572 // Чтобы увидеть решение f для v=40, disp(sol(1,:)).
2573 // Обратите внимание на то, что оно исчезает точно один раз в
2574 // [0,1] в x близком 0.98, и становится очень малым, когда
2575 // x -> 0 и очень большим, когда x -> y.
2576 // Это явно отличается от случая при малом v.
2577 // Процедура продолжения позволяет исследовать эти экспоненциальные
2578 // поведения без перехода в другие собственные состояния.
2580 ]]></programlisting>
2584 <refsection role="see also">
2585 <title>Смотрите также</title>
2586 <simplelist type="inline">
2588 <link linkend="link">link</link>
2591 <link linkend="external">external</link>
2594 <link linkend="ode">ode</link>
2597 <link linkend="dassl">dassl</link>
2602 <title>Используемые функции</title>
2603 <para>Эта функция основана на процедуре Fortran
2604 <literal>colnew</literal>, разработанной
2606 <para>U. Ascher, Department of Computer Science, University of British
2607 Columbia, Vancouver, B.C. V6T 1W5, Canada
2609 <para>G. Bader, institut f. Angewandte mathematik university of
2610 Heidelberg; im Neuenheimer feld 294d-6900 Heidelberg 1
2614 <title>Литература</title>
2617 <para>U. Ascher, J. Christiansen and R.D. Russell, collocation
2618 software for boundary-value ODEs, acm trans. math software 7 (1981),
2619 209-222. this paper contains EXAMPLES where use of the code is
2624 <para>G. Bader and U. Ascher, a new basis implementation for a mixed
2625 order boundary value ode solver, siam j. scient. stat. comput.
2630 <para>U. Ascher, J. Christiansen and R.D. Russell, a collocation
2631 solver for mixed order systems of boundary value problems, math. comp.
2636 <para>U. Ascher, J. Christiansen and R.D. russell, colsys - a
2637 collocation code for boundary value problems, lecture notes comp.sc.
2638 76, springer verlag, b. children et. al. (eds.) (1979), 164-185.
2642 <para>C. Deboor and R. Weiss, solveblok: a package for solving almost
2643 block diagonal linear systems, acm trans. math. software 6 (1980),