Получение аналитического приближенного решения обыкновенного дифференциального уравнения в виде частичной суммы ряда Тейлора функции решения

Получение аналитического приближенного решения обыкновенного дифференциального уравнения в виде частичной суммы ряда Тейлора функции решения Статьи

Рассмотрим некоторое обыкновенное дифференциальное уравнение (ОДУ), скажем \(y'+y=x\). С помощью встроенной функции DSolve можно с легкостью найти его общее аналитическое решение:

In[1]:=
DSolve[y'[x]+y[x]==x, y, x]
Out[1]=
Получение аналитического приближенного решения обыкновенного дифференциального уравнения в виде частичной суммы ряда Тейлора функции решения

или решение задачи Коши:

In[2]:=
exactSolution=DSolve[{y'[x]+y[x]==x, y[0]==1}, y, x]
Out[2]=
Получение аналитического приближенного решения обыкновенного дифференциального уравнения в виде частичной суммы ряда Тейлора функции решения

А с помощью функции NDSolve — численное решение задачи Коши:

In[3]:=
numericalSolution=NDSolve[{y'[x]+y[x]==x, y[0]==1}, y, {x, 0, 2}]
Out[3]=
Получение аналитического приближенного решения обыкновенного дифференциального уравнения в виде частичной суммы ряда Тейлора функции решения

Можно визуализировать полученные решения — аналитическое и численное, а также сравнить их:

In[4]:=
Panel@Grid[{{Plot[Evaluate[{y[x]/.exactSolution, y[x]/.numericalSolution}], {x, 0, 2}, ImageSize->500, PlotStyle->Thick, PlotLabel->Style["yаналитич. и yчисл.", 20, FontFamily->"Myriad Pro Cond", Bold], TicksStyle->Directive[Bold, 14, FontFamily->"Myriad Pro Cond"], AxesStyle->Arrowheads[{0, 0.03}], AxesLabel->(Style[#, Bold, 18, Italic, FontFamily->"Myriad Pro Cond"]&/@{"x", "y"})]}, {Plot[-Log10@Abs@Evaluate[(y[x]/.exactSolution)-(y[x]/.numericalSolution)], {x, 0, 2}, ImageSize->500, PlotPoints->1000, PlotStyle->Thick, PlotLabel->Style["-lg|yаналитич.-yчисл.|", 20, FontFamily->"Myriad Pro Cond", Bold], TicksStyle->Directive[Bold, 14, FontFamily->"Myriad Pro Cond"], AxesStyle->Arrowheads[{0, 0.03}], AxesLabel->(Style[#, Bold, 18, Italic, FontFamily->"Myriad Pro Cond"]&/@{"x", "y"})]}}, Frame->All]
Out[4]=
Получение аналитического приближенного решения обыкновенного дифференциального уравнения в виде частичной суммы ряда Тейлора функции решения

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

Для него нам известно начальное условие \(y(0)=1\). Подставив значение \(x=0\) в уравнение, получим:

\(y'(0)+y(0)=0\), откуда следует, что \(y'(0)=-y(0)\), т. е. \(y'(0)=-1\).

Продифференцировав ОДУ по переменной \(x\), получим:

\(y\text{''}(x)+y'(x)=1\).

Подставив в это уравнение \(x=0\), получим:

\(y\text{''}(0)+y'(0)=1\), или \(y\text{''}(0)=1-y'(0)\), т. е. \(y\text{''}(0)=2\).

Продолжая этот процесс, получим, что:

\(y\text{'''}(0)=-2\), \(y\text{''''}(0)=2\), \(y^{(5)}(0)=-2\), \(y^{(6)}(0)=2\),

\(y^{(7)}(0)=-2\), \(y^{(8)}(0)=2\), ... .

Это означает, что решение может быть записано в виде:

\(y(x)=\underset{k=0}{\overset{+\infty}{\Sigma}}\frac{x^k}{k!}y^{(k)}(0)=\frac{x^0}{0!}1+\frac{x^1}{1!}(-1)+\underset{k=2}{\overset{+\infty}{\Sigma}}\frac{x^k}{k!}2(-1)^k=\)

\(=1-x+2\underset{k=2}{\overset{+\infty}{\Sigma}}\frac{(-1)^kx^k}{k!}.\)

Полученный ряд является решением ОДУ, а также может служить аналитическим приближением к решению. Посмотрим, как ведет себя полученное аналитическое приближение:

In[5]:=
Manipulate[Grid[{{Plot[Evaluate[{E-x (2-Ex+Ex x), 1-x+2Sum[(-1)kxkk!, {k, 2, kol}]}], {x, 0, 5}, PlotRange->{0, 5}, PlotStyle->Thick, ImageSize->500, AspectRatio->1/2, PlotLegends->Placed[{Style["yаналитич.", 20, FontFamily->"Myriad Pro Cond", Bold], Style["yаналитич. прибл.", 20, FontFamily->"Myriad Pro Cond", Bold]}, Above], TicksStyle->Directive[Bold, 14, FontFamily->"Myriad Pro Cond"], AxesStyle->Arrowheads[{0, 0.03}], AxesLabel->(Style[#, Bold, 18, Italic, FontFamily->"Myriad Pro Cond"]&/@{"x", "y"})]}, {Plot[Evaluate[-Log10@Abs@{E-x (2-Ex+Ex x)-(1-x+2Sum[(-1)kxkk!, {k, 2, kol}])}], {x, 0, 5}, PlotStyle->Thick, AspectRatio->1/2, ImageSize->500, PlotRange->{-1, 20}, WorkingPrecision->30, PlotLabel->Style["-lg|yаналитич.-yаналитич. прибл.|", 20, FontFamily->"Myriad Pro Cond", Bold], TicksStyle->Directive[Bold, 14, FontFamily->"Myriad Pro Cond"], AxesStyle->Arrowheads[{0, 0.03}], AxesLabel->(Style[#, Bold, 18, Italic, FontFamily->"Myriad Pro Cond"]&/@{"x", "y"})]}}, Frame->All], {{kol, 2, "Максимальное значение k:"}, 2, 20, 1, Appearance->"Open"}]
Out[5]=
Получение аналитического приближенного решения обыкновенного дифференциального уравнения в виде частичной суммы ряда Тейлора функции решения

В Mathematica вы можете найти сумму ряда в правой части полученного аналитического приближения:

In[6]:=
1-x+2Sum[(-1)kxkk!, {k, 2, Infinity}]
Out[6]=
Получение аналитического приближенного решения обыкновенного дифференциального уравнения в виде частичной суммы ряда Тейлора функции решения

И убедиться, что решение в виде ряда Тейлора (Маклорена в данном случае), сходится к точному решению:

In[7]:=
{Simplify[%], Simplify[y[x]/.exactSolution][[1]]}
Out[7]=
Получение аналитического приближенного решения обыкновенного дифференциального уравнения в виде частичной суммы ряда Тейлора функции решения

Теперь создадим функцию, которая проделывает все описанное выше:

In[8]:=
ODESeriesSolve[{ode_, initials__}, fName_, {variableName_, n_}, OptionsPattern[FindExactSolution->False]]:=Block[{system, derivatives, whatToFind, order, initialX, seriesTerms}, order=Max[Cases[FullForm[ode], Derivative[i_][_][_]:>i, Infinity]];

initialX=Cases[{initials}, fName[x_]:>x, Infinity][[1]];

system=Table[D[ode, {x, i}]/.variableName->initialX, {i, 0, n}];

whatToFind=DeleteDuplicates[Cases[FullForm[system], Derivative[_][_][_], Infinity]~Join~Cases[{initials}, Equal[x_, _]:>x, Infinity]];

derivatives=Solve[{initials}~Join~system, whatToFind][[1]];

seriesTerms=Table[(D[fName[variableName], {variableName, i}]/.variableName->initialX)*(variableName-initialX)^i/i!, {i, 0, n+order}]/.derivatives;

If[OptionValue[FindExactSolution], {Sort[derivatives], seriesTerms, DSolve[{ode, initials}, fName, variableName]}, {Sort[derivatives], seriesTerms}]]

Можем применить ее к решению разнообразных ОДУ и посмотреть что она предоставляет нам, как в частном:

In[9]:=
ODESeriesSolve[{y''[x]+x y[x]==0, y[0]==1, y'[0]==0}, y, {x, 20}, FindExactSolution->True]
Out[9]=
Получение аналитического приближенного решения обыкновенного дифференциального уравнения в виде частичной суммы ряда Тейлора функции решения

так и в общем виде:

In[10]:=
ODESeriesSolve[{y''[x]+x y[x]==A, y[x0]==y0, y'[x0]==y1}, y, {x, 20}, FindExactSolution->True]
Out[10]=
Получение аналитического приближенного решения обыкновенного дифференциального уравнения в виде частичной суммы ряда Тейлора функции решения

Если есть желание, можно посмотреть на то, как последовательные аналитические приближения описывают решение:

In[11]:=
{coeff, terms, solution}=ODESeriesSolve[{y''[x]+x y[x]==0, y[0]==1, y'[0]==0}, y, {x, 20}, FindExactSolution->True]
Out[11]=
Получение аналитического приближенного решения обыкновенного дифференциального уравнения в виде частичной суммы ряда Тейлора функции решения
In[12]:=
curves=DeleteDuplicates[Accumulate[terms]];

exact=y[x]/.solution;

Plot[{curves, exact}, {x, 0, 6}, PlotRange->{-1, 2}, ImageSize->600, AspectRatio->Automatic, PlotStyle->Table[{Thick, RGBColor[i/Length[curves], 0, 0]}, {i, 1, Length[curves]}]~Join~{{Blue, Thick}}, TicksStyle->Directive[Bold, 14, FontFamily->"Myriad Pro Cond"], AxesStyle->Arrowheads[{0, 0.03}], AxesLabel->(Style[#, Bold, 18, Italic, FontFamily->"Myriad Pro Cond"]&/@{"x", "y"})]
Out[12]=
Получение аналитического приближенного решения обыкновенного дифференциального уравнения в виде частичной суммы ряда Тейлора функции решения

В тех случаях, когда получить аналитическое решение не представляется возможным, такой метод позволит получить аналитическое приближение:

In[13]:=
{coeff, terms}=ODESeriesSolve[{y'[x]==y[x]^2+ x-Cos[2x], y[0]==0}, y, {x, 20}, FindExactSolution->False]
Out[13]=
Получение аналитического приближенного решения обыкновенного дифференциального уравнения в виде частичной суммы ряда Тейлора функции решения
In[14]:=
curves=DeleteDuplicates[Accumulate[terms]];

numerical=y[x]/.NDSolve[{y'[x]==y[x]^2+ x-Cos[2x] , y[0]==0}, y, {x, 0, 1.5}];

Plot[{curves, numerical}, {x, 0, 1.5}, PlotRange->{-5, 5}, ImageSize->600, AspectRatio->1, PlotStyle->Table[{Thick, RGBColor[i/Length[curves], 0, 0]}, {i, 1, Length[curves]}]~Join~{{Thick, Blue}}, TicksStyle->Directive[Bold, 14, FontFamily->"Myriad Pro Cond"], AxesStyle->Arrowheads[{0, 0.03}], AxesLabel->(Style[#, Bold, 18, Italic, FontFamily->"Myriad Pro Cond"]&/@{"x", "y"})]
Out[14]=
Получение аналитического приближенного решения обыкновенного дифференциального уравнения в виде частичной суммы ряда Тейлора функции решения

С помощью функции Manipulate можно продемонстрировать приближение к решению с помощью частичных сумм ряда Тейлора в интерактивном виде:

In[15]:=
Manipulate[Plot[{curves$, curves$[[n]], numerical$}, {x, 0, 1.5}, PlotRange->{-5, 5}, ImageSize->500, AspectRatio->1, PlotStyle->Table[Gray, {i, 1, Length[curves$]}]~Join~{{AbsoluteThickness[5], Orange}, {AbsoluteThickness[3], Blue}}, TicksStyle->Directive[Bold, 14, FontFamily->"Myriad Pro Cond"], AxesStyle->Arrowheads[{0, 0.03}], AxesLabel->(Style[#, Bold, 18, Italic, FontFamily->"Myriad Pro Cond"]&/@{"x", "y"})], {{n, 1, "Количество членов ряда:"}, 1, Length[curves$], 1}, Initialization:>{terms$={0, -x, x22, x3, -x44, -29 x560, x64, 65 x7252, -89 x8480, -847 x96480, 6907 x1050400, 317869 x114989600, -37139 x12388800, -1479463 x1355598400, 5403121 x1483825280, 318103439 x1540864824000, -832036307 x1619813248000, 11323481131 x1711115232128000, 7431184181 x18280215936000, -185847716111 x1942237882086400, -327637807877 x2020209512960000, 2031473513846531 x21399147985716480000};

curves$=DeleteDuplicates[Accumulate[terms$]];

numerical$=y[x]/.NDSolve[{y'[x]==y[x]^2+ x-Cos[2x] , y[0]==0}, y, {x, 0, 1.5}];

}]
Out[15]=
Получение аналитического приближенного решения обыкновенного дифференциального уравнения в виде частичной суммы ряда Тейлора функции решения

Оценить статью
Блог о Wolfram Mathematica

Оставить комментарий

avatar
  Подписаться  
Уведомление о