Mathematica в действии, блог Осипова Романа

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

Нажмите, чтобы скачать данный пост в виде ноутбука Mathematica...

Общее количество использованных в посте встроенных функций или символов: 82

Список имен используемых встроенных функций и символов в порядке их появления в коде:
DSolve | Equal (==) | Plus (+) | Derivative (') | Set (=) | List ({...}) | NDSolve | Panel | Grid | Plot | Evaluate | ReplaceAll (/.) | Rule (->, ->) | ImageSize | PlotStyle | Thick | PlotLabel | Style | FontFamily | Bold | TicksStyle | Directive | AxesStyle | Arrowheads | AxesLabel | Map (/@) | Function (&) | Slot (#) | Italic | Times (*, ×) | Log10 | Abs | PlotPoints | Frame | All | Manipulate | Power (^) | E (e) | Sum (∑) | Factorial (!) | PlotRange | AspectRatio | PlotLegends | Placed | Above | WorkingPrecision | Appearance | Infinity | Simplify | Out (%) | Part ([[…]]) | SetDelayed (:=) | Pattern (:) | Blank (_) | BlankSequence (__) | OptionsPattern | False | Block | CompoundExpression (;) | Max | Cases | FullForm | RuleDelayed (:>, :->) | Table | D (∂) | DeleteDuplicates | Join | Solve | If | OptionValue | Sort | True | Accumulate | Automatic | RGBColor | Length | Blue | Cos | Gray | AbsoluteThickness | Orange | Initialization

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

закрыть

DSolve[y'[x] + y[x] == x, y, x]

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

закрыть

exactSolution = DSolve[{y'[x] + y[x] == x, y[0] == 1}, y, x]

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

закрыть

numericalSolution = NDSolve[{y'[x] + y[x] == x, y[0] == 1}, y, {x, 0, 2}]

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

закрыть

Panel@Grid[{{Plot[ Evaluate[{y[x] /. exactSolution, y[x] /. numericalSolution}], {x, 0, 2}, ImageSize -> 500, PlotStyle -> Thick, PlotLabel -> Style["\!\(\*SubscriptBox[\(y\), \(\(аналитич\)\(.\)\)]\) и \\!\(\*SubscriptBox[\(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|\!\(\*SubscriptBox[\(y\), \\(\(аналитич\)\(.\)\)]\)-\!\(\*SubscriptBox[\(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]

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

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

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

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

y''(x)+y'(x)=1.

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

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

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

y'''(0)=-2, y''''(0)=2, , ,

, , ... .

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

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

закрыть

Manipulate[ Grid[{{Plot[ Evaluate[{E^-x (2 - E^x + E^x x), 1 - x + 2 Sum[((-1)^k x^k)/k!, {k, 2, kol}]}], {x, 0, 5}, PlotRange -> {0, 5}, PlotStyle -> Thick, ImageSize -> 500, AspectRatio -> 1/2, PlotLegends -> Placed[{Style[ "\!\(\*SubscriptBox[\(y\), \(\(аналитич\)\(.\)\)]\)", 20, FontFamily -> "Myriad Pro Cond", Bold], Style["\!\(\*SubscriptBox[\(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 - E^x + E^x x) - (1 - x + 2 Sum[((-1)^k x^k)/k!, {k, 2, kol}])}], {x, 0, 5}, PlotStyle -> Thick, AspectRatio -> 1/2, ImageSize -> 500, PlotRange -> {-1, 20}, WorkingPrecision -> 30, PlotLabel -> Style["-lg|\!\(\*SubscriptBox[\(y\), \\(\(аналитич\)\(.\)\)]\)-\!\(\*SubscriptBox[\(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"}]

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

закрыть

1 - x + 2 Sum[((-1)^k x^k)/k!, {k, 2, Infinity}]

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

закрыть

{Simplify[%], Simplify[y[x] /. exactSolution][[1]]}

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

закрыть

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}]]

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

закрыть

ODESeriesSolve[{y''[x] + x y[x] == 0, y[0] == 1, y'[0] == 0}, y, {x, 20}, FindExactSolution -> True]

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

закрыть

ODESeriesSolve[{y''[x] + x y[x] == A, y[x0] == y0, y'[x0] == y1}, y, {x, 20}, FindExactSolution -> True]

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

закрыть

{coeff, terms, solution} = ODESeriesSolve[{y''[x] + x y[x] == 0, y[0] == 1, y'[0] == 0}, y, {x, 20}, FindExactSolution -> True]

закрыть

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"})]

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

закрыть

{coeff, terms} = ODESeriesSolve[{y'[x] == y[x]^2 + x - Cos[2 x], y[0] == 0}, y, {x, 20}, FindExactSolution -> False]

закрыть

curves = DeleteDuplicates[Accumulate[terms]]; numerical = y[x] /. NDSolve[{y'[x] == y[x]^2 + x - Cos[2 x] , 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"})]

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

закрыть

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, x^2/2, x^3, -(x^4/4), -((29 x^5)/60), x^6/4, (65 x^7)/ 252, -((89 x^8)/480), -((847 x^9)/6480), (6907 x^10)/50400, ( 317869 x^11)/ 4989600, -((37139 x^12)/388800), -((1479463 x^13)/55598400), ( 5403121 x^14)/83825280, (318103439 x^15)/ 40864824000, -((832036307 x^16)/19813248000), (11323481131 x^17)/ 11115232128000, (7431184181 x^18)/ 280215936000, -((185847716111 x^19)/42237882086400), -(( 327637807877 x^20)/20209512960000), (2031473513846531 x^21)/ 399147985716480000}; curves$ = DeleteDuplicates[Accumulate[terms$]]; numerical$ = y[x] /. NDSolve[{y'[x] == y[x]^2 + x - Cos[2 x] , y[0] == 0}, y, {x, 0, 1.5}];}]

Блог принадлежит “Русскоязычной поддержке Wolfram Mathematica"©
При любом использовании материалов блога, ссылка на блог обязательна.
Создано с помощью Wolfram Mathematica 9