NDSolve podaje złe rozwiązanie

Question

Rozważmy ODE $$\frac{y'y}{1+\frac{1}{2} \sqrt{1+ y'^2}}=-x.$$

Za pomocą

NDSolve[{-x==y'[x] y[x]/(1+Sqrt[1+(y'[x])^2]/2), y[0] ==3}, y, {x,-7,7}]

i kreślenie prowadzi do dwóch rozwiązań (niektóre ostrzeżenia w pobliżu granic)

$$y(x) = \sqrt{6^2 - x^2} - 3$$ i $$y(x) = \sqrt{2^2 - x^2} + 1.$$

Ale tylko to drugie jest poprawnym rozwiązaniem! Bez względu na to, którą „metodę” próbowałem, zawsze otrzymywałem całkowicie błędne rozwiązanie. Z wyjątkiem używania

Method -> {"EquationSimplification" -> "Residual"}

Dlaczego?

Uwaga: Jak wskazano w odpowiedzi poniżej, ustalenie wartości w x=0jest krytyczne, ponieważ$y'$znika tutaj. Ale użycie innych wartości początkowych, takich jak y[Sqrt[3]]=2problem, staje się jeszcze gorsze, ponieważ jedna gałąź jest teraz całkowicie błędna wszędzie, a druga gałąź jest poprawna tylko na małym obszarze.

Answer 1

Przyczyną tego zachowania wydaje się być duży błąd logiczny w programie NDSolve. Wydaje się, że podczas obliczeń traktuje wyrażenia takie jak: y==Sqrt[x]i y^2==xtak samo. Ale, jak każdy użytkownik tutaj wie, tak nie jest!

Jako potwierdzenie weź swój konkretny przykład: Mnożenie przez mianownik daje $$-x\left(1-\frac{1}{2} \sqrt{1+(y'(x))^2}\right)=y'(x) y(x).$$ Głupie rozwiązywanie obu stron i rozwiązywanie problemów $y'(x)$ tworzy dwie gałęzie

NDSolve[{y'[x]==(4 x y[x]+Sqrt[3 x^4 + 4 x^2 y[x]^2])/(x^2 - 4 y[x]^2) , y[0]==3}, y, {x,-6,6}]

i

NDSolve[{y'[x]==(4 x y[x]-Sqrt[3 x^4 + 4 x^2 y[x]^2])/(x^2 - 4 y[x]^2) , y[0]==3}, y, {x,-6,6}]

Są to rzeczywiście dokładnie te gałęzie NDSolve, które zapewniają, chociaż żadna nie jest ważna.

Co gorsza, choć jest fundamentalny, nie sprawdza rozwiązań. Wymagałoby to tylko dodatkowej linii kodu w algorytmie, ponieważ używa on już krotek$(x_i,y(x_i),y'(x_i)$. Po prostu podłącz je do równania i sprawdź, czy jest prawdziwe, czy fałszywe (do pewnego błędu numerycznego).

Edycja: NDSolve musi przekształcić równanie do jakiejś standardowej postaci, którą kontroluje EquationSimplification. Istnieją trzy możliwe warianty tej metody: MassMatrix, Residuali Solvektóry jest domyślnym. Ten ostatni przekształca równanie w postać bez pochodnych po jednej stronie. Układ jest następnie rozwiązywany za pomocą zwykłego narzędzia do rozwiązywania równań różniczkowych. Gdy Residualzostanie wybrany, wszystkie niezerowe wyrazy w równaniu są po prostu przenoszone na jedną stronę, a następnie rozwiązywane za pomocą różniczkowego równania algebraicznego . To jest powód, dla którego wynik jest poprawny w tym przypadku, ponieważ nie używa tego, Solveco jest tutaj błędne .

Answer 2

Ogólny problem

Używając NDSolvedo rozwiązywania IVP pierwszego rzędu, istnieją zasadniczo dwa sposoby skonfigurowania ODE:

y'[x] == f[x, y[x]]     (* explicit form *)
F[x, y[x], y'[x]] == 0  (* implicit form *)

Większość rozwiązywania numerycznych wymaga określenia problemu w jawnej formie. W Mathematica jest tylko jeden solwer, który działa z formą niejawną, IDA , i jest ograniczony do precyzji maszyny. Ponieważ łatwo jest przekształcić niejawną postać w jawną postać równania drugiego rzędu poprzez różnicowanie względem x, być może nie było dużego nacisku na opracowanie solwerów z formą niejawną.

W Mathematica możesz poprosić o wypróbowanie rozwiązania w dowolnej formie z Methodopcją:

Method -> {"EquationSimplification" -> "Solve"}    (* explicit *)
Method -> {"EquationSimplification" -> "Residual"} (* implicit *)

Z "Solve"metodą, która jest domyślna, NDSolvewywołuje Solvekonwersję ODE do postaci jawnej. Równanie podane w postaci niejawnej może mieć wiele rozwiązań, a jeśli tak, NDSolvebędzie całkować każde z nich osobno. Tak dzieje się w przykładzie PO. Ponadto NDSolvejest skonfigurowany do niezależnej integracji oddzielnych jawnych form ODE i nie może ich łączyć, co jest wymagane w przypadku OP (patrz odpowiedź @ BobHanlon ).

Obecnie Solvekwestia hojności odgrywa tutaj ważną rolę. W przypadku OP zwraca rozwiązania, które są ważne w niektórych domenach i nieważne w innych niepustych regionach, w tym w tych, które chcemy zintegrować. Reducejest znacznie ostrożniejszy i poprawnie analizuje system PO. Można Solveskorzystać Reducez tej opcji Method -> Reduce, ale nadal zwraca ona dwa oddzielne rozwiązania, każde poprawne po jednej stronie x == 0. Ponadto zwraca ConditionalExpression, który NDSolvedławi na (i daje „non-numerycznego” NDSolve::ndnumbłąd w początkowym stanie podczas tej ProcessEquationsfazy ). ConditionalExpressionzostał wprowadzony dość późno, w wersji 8, i być może nie dostało wystarczającej liczby żądań, aby NDSolveprawidłowo go obsłużyć, do WRI.

OTOH, "Residual"metoda rozwiązuje ODE niejawnie na każdym kroku. Ponieważ oba rozwiązania są ważne jednocześnie tylko pod adresem x == 0, po NDSolvewykonaniu kroku zostanie znaleziona właściwa gałąź . To oblicza poprawne rozwiązanie, o którym wspomina PO. Jedyną wadą jest to, że dostępna jest tylko jedna metoda integracji i tylko w precyzji maszyny.

Wydaje się, że na NDSolve`ProcessEquationsetapie łatwo byłoby sprawdzić, czy pierwotna niejawna forma ODE jest spełniana przez formy jawne w warunku początkowym. To nie uchwyciłoby problemu z przykładu w y[0] == 3punkcie, w którym obie formy jawne spełniają niejawną postać ODE, ale wychwyciłyby problem w y[1] == 2. Inną kwestią związaną z rozwiązaniami zwracanymi przez Solvejest to, że jawna formuła na y'[x]potrzeby przełączania gałęzi na inne rozwiązanie zwracana jest Solvepo przecięciu integracji x == 0. Przełączanie gałęzi nie jest czymś, co NDSolvejest przewidziane, ani nie wydaje mi się łatwą poprawką programistyczną, ponieważ każde rozwiązanie jest zintegrowane niezależnie. Poniżej przedstawiono kilka sposobów, aby to zrobić, ale wszystkie wymagają od użytkownika przygotowania NDSolvepołączenia. Żadne nie są wykonywane automatycznie przez NDSolve, co byłoby pożądane.

Wreszcie, czego powinien oczekiwać użytkownik? W obliczeniach naukowych od dawna oczekiwano, że użytkownik ustawi numeryczną całkowanie równań różniczkowych. Wydaje się, że nadal tak jest w MATLAB i NumPy. Nie znam klonu wystarczająco dobrze, żeby to skomentować. Ogólna filozofia Mathematica polega na tym, aby wszystko było jak najbardziej zautomatyzowane. Mathematica miała również tendencję do stosowania ogólnie rzecz biorąc prawdziwych rozwiązań zamiast bardziej rygorystycznych ograniczeń. Są one tu nieco sprzeczne, ponieważ ogólne metody Solvesą źródłem problemu z NDSolverozwiązaniami. Z drugiej strony, aby wszystko odbywało się automatycznie, jest nie tyle celem Wolframa, co zasadą przewodnią. Pytania i odpowiedzi na tej stronie pokazują, że Automaticnie zawsze spełnia to swoje zadanie. Użytkownik często musi zrozumieć problem, wiedzieć, jakie rozwiązania są dostępne, odpowiednio przygotować dane wejściowe i zadzwonić do solvera z odpowiednimi opcjami. W przypadku niejawnego IVP użytkownik powinien mieć świadomość, że może wystąpić problem z rozwiązaniem dla y'[x]. Powinni również być świadomi, że istnieją standardowe sposoby radzenia sobie z niejawnymi formami ODE:

• użycie niejawnego solvera, takiego jak IDA, wywoływanego, gdy "Residual"jest wywoływany;
• różnicowanie w celu podniesienia zamówienia;
• rozwiązanie y'[x]jawnie, "Solve"metoda domyślna .

Powtórzę, że moim zdaniem rozsądne jest oczekiwanie NDSolvesprawdzenia, czy jawna forma spełnia pierwotną niejawną postać ODE w początkowym stanie. Chociaż użytkownik może sprawdzić wyniki NDSolvepo fakcie, w przypadkach takich jak IVP y[1] == 2zapobiegnie to obcej integracji.

Przykłady PO

Jawne rozwiązania dla y'[x]ODE PO mają dwie gałęzie dla x < 0i dwie dla x > 0. Oba rozwiązania wynikają z (algebraicznej) racjonalizacji zmiennej różniczkowej, która wprowadza możliwość rozwiązań obcych. W rzeczywistości zestaw rozwiązań składa się z czterech połączonych komponentów, dwóch dla interwału x < 0i dwóch dla x > 0. Każde rozwiązanie zwrócone przez Solvejest ważne w jednym przedziale, ale nie w obu. Możemy jednak przekształcić je w jedno poprawne i jedno nieprawidłowe rozwiązanie Simplify[..., x > 0], ale myślę, że nie jest to ogólna technika.

Obejście nr 1

Odkrycie PO:

ode = -x == y'[x] y[x]/(1 + Sqrt[1 + (y'[x])^2]/2);

ListLinePlot[
 NDSolveValue[{ode, y[0] == 3}, y, {x, -7, 7}, 
  Method -> {"EquationSimplification" -> "Residual"}],
 PlotRange -> All
 ]

Obejście nr 2

Różnicowanie ODE podnosi kolejność, ale skutkuje taką, dla której istnieje unikalna jawna forma. Musisz użyć ODE, aby rozwiązać warunek początkowy dla y'[0].

sol = NDSolve[{D[ode, x], y[0] == 3, y'[0] == 0}, y, {x, -7, 7}]

Obejście nr 3

Użyj poprawnej jawnej formy , zbudowanej z właściwych gałęzi dla x <> 0:

ode2 = y'[x] == 
   Piecewise[{
    {(4 x y[x] - Sqrt[3 x^4 + 4 x^2 y[x]^2])/(x^2 - 4 y[x]^2), x < 0}}, 
    (4 x y[x] + Sqrt[3 x^4 + 4 x^2 y[x]^2])/(x^2 - 4 y[x]^2)];

sol = NDSolve[{ode2, y[0] == 3}, y, {x, -7, 7}]

Obejście nr 4

Istnieją problemy z naszą notacją algebraiczną i jej związkiem z funkcjami algebraicznymi. Zastosowanie tego założenia x > 0 zmienia wybór cięcia gałęzi przy upraszczaniu zwracanych rozwiązań, Solvetak aby jedno z nich było poprawne. Innymi słowy, daje to prostszą formułę, y'[x]która jest odpowiednikiem Obejścia nr 3.

sol = NDSolve[{#, y[0] == 3} /. Rule -> Equal, y, {x, -7, 7}] & /@
  Assuming[x > 0,
   Select[Simplify@Solve[ode, y'[x]], 
    ode /. # /. {y[x] -> 1, x -> 1.`20} &]
   ] // Apply[Join]

Obejście nr 5

Ta Solveopcja Method -> Reducetworzy poprawne rozwiązania w postaci ConditionalExpression. Aby uzyskać metodę, która sprawdza i wybiera poprawną gałąź ODE, która niejawnie definiuje y'[x], użytkownik musiałby wykonać własne przetwarzanie wstępne. Poniżej przedstawiono sposób, w jaki rhs[]wybiera się gałąź, która spełnia oryginalne ODE, poprzez konwersję wyrażeń warunkowych na pojedynczą Piecewisefunkcję. Warunki są konwertowane z równań a == bna porównanie Abs[a-b] < 10^-8. Musiałem x == 0ręcznie dodać wartość w punkcie rozgałęzienia .

Innymi słowy, sprawdza to y'[x]na każdym kroku i wybiera właściwą gałąź dla tego kroku. W ten sposób automatycznie przełącza gałęzie w razie potrzeby, x == 0w przypadku problemu PO. Warto zwrócić uwagę, że rozwiązuje to problem wynikający z racjonalizacji ODE wprowadzającej obce gałęzie. Możliwe jest, że niejawny formularz ODE ma wiele prawidłowych gałęzi. Poniższa metoda połączy je wszystkie (jeśli rozwiązania mają ConditionalExpressionformę), co należy uznać za błąd, chociaż nadal może przypadkowo dać prawidłowe rozwiązanie. Dla ODE PO robi to dobrze.

ClearAll[rhs];
rhs[x_?NumericQ, y_?NumericQ] = Piecewise[
   yp /. Solve[ode /. {y[x] -> y, y'[x] -> yp}, yp, 
       Method -> Reduce] /. ConditionalExpression -> List /. 
    Equal -> (Abs[#1 - #2] < 10^-8 &),
   0 (* y'[0] == 0 *)];

sol = NDSolve[{y'[x] == rhs[x, y[x]], y[0] == 3}, y, {x, -7, 7}]

Oto bardzo hakerski sposób na naprawienie wyniku wewnętrznego Solve. Osiąga się to za pomocą sekwencji wirusowych UpValuesdo $tagże przepisuje się ConditionalExpressionroztwór do Piecewiseroztworu, jak to powyżej.

opts = Options@Solve;
SetOptions[Solve, Method -> Reduce];
Block[{ConditionalExpression = $tag, $tag},
 $tag /: Rule[v_, $tag[a_, b_]] := $tag[v, a, b];
 $tag /: {$tag[v_, a_, b_]} := $tag[List, v, a, b];
 $tag /: call : {$tag[List, v_, __] ..} := {{v -> 
     Piecewise[
      Unevaluated[call][[All, -2 ;;]] /. $tag -> List /. 
       Equal -> (Abs[#1 - #2] < 1*^-8 &)]}};
 sol = NDSolve[{ode, y[0] == 3}, y, {x, -7, 7}]
 ]
SetOptions[Solve, opts];

Jak zobaczyć, co Solverobi w środkuNDSolve

Jeśli chcesz zobaczyć, co dzieje się wewnętrznie, możesz użyć Trace. NDSolveużywa Solverozwiązania ODE dla pochodnej najwyższego rzędu, jeśli to możliwe, i używa rozwiązania (rozwiązań) do skonstruowania całki (całek). To pokazuje Solvewywołanie i jego wartość zwracaną:

Trace[
 NDSolve[
  {ode, y[0] == 3},
  y, {x, -7, 7}],
 _Solve,
 TraceForward -> True,
 TraceInternal -> True
 ]

Answer 3

Clear["Global`*"]

sol = DSolve[{-x == y'[x] y[x]/(1 + Sqrt[1 + (y'[x])^2]/2), y[0] == 3}, y, 
   x] // Quiet

(* {{y -> Function[{x}, Sqrt[5 - x^2 + 2 Sqrt[4 - x^2]]]}, 
    {y -> Function[{x}, Sqrt[45 - x^2 - 6 Sqrt[36 - x^2]]]}} *)

FunctionDomain[y[x] /. sol[[1]], x]

(* -2 <= x <= 2 *)

Pierwsze rozwiązanie dotyczy -2 <= x <= 2

{-x == y'[x] y[x]/(1 + Sqrt[1 + (y'[x])^2]/2), y[0] == 3} /. sol[[1]] // 
 Simplify[#, -2 <= x <= 2] &

(* {True, True} *)

FunctionDomain[y[x] /. sol[[2]], x]

(* -6 <= x <= 6 *)

Drugie rozwiązanie jest prawdziwe dla x == 0

{-x == y'[x] y[x]/(1 + Sqrt[1 + (y'[x])^2]/2), y[0] == 3} /. sol[[2]] // 
 FullSimplify[#, -6 <= x <= 6] &

(* {x == 0, True} *)

Plot[Evaluate[y[x] /. sol], {x, -6, 6},
 PlotLegends -> Placed[Automatic, {.75, .2}]]

W przypadku rozwiązania numerycznego ogranicz domenę do {- 2, 2}

soln = NDSolve[{-x == y'[x] y[x]/(1 + Sqrt[1 + (y'[x])^2]/2), y[0] == 3}, 
    y, {x, -2, 2}] // Quiet;

Rozwiązania numeryczne obowiązują w różnych częściach domeny

Plot[Evaluate[y[x] /. soln], {x, -2, 2},
 PlotRange -> {0, 3.1},
 PlotLegends -> Placed[Automatic, {.7, .5}]]

Answer 4

Zacznij od

Plot[Evaluate[y[x] /. sol], {x, -2, 2}, 
 PlotLegends -> Placed[Automatic, {.75, .2}], PlotPoints -> 1600, 
 ImageSize -> Large, PlotRange -> Full]

Co jest w równaniu różniczkowym?

$$\frac{𝑦′𝑦}{1+\sqrt{1+𝑦′^2}}=−𝑥$$

1. To równanie różniczkowe typu niejawnego.

2. Jest to równanie różniczkowe pierwszego rzędu ${y,y'}$.

3. Jest to nieliniowe równanie różniczkowe.

4. Podaje się go w postaci ilorazu, więc trzeba szukać osobliwości mianownika.

5. W mianowniku, który ma być traktowany, dokonuje się wyboru znaku pierwiastka drugiego stopnia. W rzeczywistości mianownik nie może wynosić zero$x$ i $y'$ tak długo, jak wybrany jest znak korzenia.

6. Istnieje postać podanego równania różniczkowego gdzie $f(x,y,y')==0$:

   y '[x] == Piecewise [{{(4 xy [x] - Sqrt [3 x ^ 4 + 4 x ^ 2 y [x] ^ 2]) / (x ^ 2 - 4 y [x] ^ 2 ), x <0}}, (4 xy [x] + Sqrt [3 x ^ 4 + 4 x ^ 2 y [x] ^ 2]) / (x ^ 2 - 4 y [x] ^ 2)]

Dzięki temu wiemy różne fakty o tym, co Mathematica może dla nas zrobić!

A. Rozwiązanie jest możliwe dzięki DSolve! DSolve rozwiązuje równanie różniczkowe dla funkcji u ze zmienną niezależną $x$ dla $x$między Subscript[x, min]a Subscript[x, max]. B. W ogóle nie potrzebujemy NDSolve. C. Ponieważ zależność funkcjonalna jest stała i różniczkowalna w danym przedziale, rozwiązanie ma te właściwości również w przedziale.

Z pytania wynika, że ​​jeden problem jest otwarty na właściwe rozwiązanie. Czym są$x_min$ i $x_max$?

Z rozwiązania DSolve:

sol = DSolve[{-x == y'[x] y[x]/(1 + Sqrt[1 + (y'[x])^2]/2), 
   y[0] == 3}, y, x]

( {{y -> Funkcja [{x}, Sqrt [5 - x ^ 2 + 2 Sqrt [4 - x ^ 2]]]}, {y -> Funkcja [{x}, Sqrt [45 - x ^ 2 - 6 Sqrt [36 - x ^ 2]]]}} )

Otrzymujemy informację, że rozwiązania nie są ograniczone do domeny rozwiązań. Mając pierwotne równanie różniczkowe jako dane wejściowe, otrzymujemy informacje, które DSolveodwołują się do wbudowanej w Mathematica metodologii obliczania odwrotnej funkcji równania różniczkowego. Dlatego przywołuje Reduce. Dane wyjściowe nie zawierają żadnych wyników z Reduce.

Są to komunikaty generowane w celu dalszego zatrzymywania takich komunikatów, jak wcześniej w pośredniej cue wyjściowej komunikatu. W końcu znajduje "obejście" # 3 z @ michael-e2, ale jest to proces wbudowany, a nie "obejście", w przeciwnym razie zestaw rozwiązań byłby pusty.

Zatem tym, co ogranicza rozwiązanie dla domeny, jest wybór pokazany przez @ bob-hanlon przy użyciu FunctionDomain. FunctionDomainogranicza się do Reals. Tego nie ma w pytaniu. A NDSolve nie ograniczyłby metod rozwiązania do Reals. Jak pokazuje moje zdjęcie wstępne, nie ma problemu z pierwszym rozwiązaniem.

Potrzebujemy pewnych rozważań geometrycznych. Podane równanie różniczkowe, nieliniowe, opisuje przesunięcia elips i tylko ich granicę. Zatem gałęzie pokazane przez @ bob-hanlon poza ograniczeniem do Realspojawiania się nie są już poprawne. Elipsy nie są nieskończenie wydłużane.

Roztwór należy poddawać dalszej obróbce, aż ocena stanie się sensowna. Matematyka wymaga, aby opisać korzenie. Nie chcemy inwersji$x(y)$. W matematyce istnieje wiele opisów elips.

Rozwiązania:

GraphicsGrid[{{Graphics[Circle[{0, 0}, {2, 3}], Axes -> True, 
    PlotRange -> {{-6, 6}, {-3.1, 3.1}}], 
   Graphics[Circle[{0, 0}, {5.2, 3}], Axes -> True, 
    PlotRange -> {{-6, 6}, {-3.1, 3.1}}]}}]

Dlaczego to mamy? OK. Wynika to z nieliniowości równania różniczkowego i samego równania różniczkowego Reals.$x(0)==3$całkowicie naprawia wielokropek. Do rozwiązania jest tylko jeden parametr. Mathematica oblicza to przy użyciu Reduce. Możemy to zrobić ręcznie, jak pokazuje inna odpowiedź. To jest konieczne.

Ten krok jest tak skomplikowany, jak zaakceptowanie, że Mathematica klasyfikuje tak, jak ja jawnie wewnętrznie przeprowadziłem równanie różniczkowe w NDSolve. Metodologia rozwiązań przekazuje proces rozwiązywania równań różniczkowych do, DSolvea następnie interpoluje rozwiązanie wzięte z tego procesu i daje wynik. Jest to szczególny przypadek oceny lenistwa. Więc moja odpowiedź nie polega na rozwiązaniu tego za pomocą, DSolveale NDSolvezamiast tego, używając ścieżki prowadzonej przez głowę.

W ten sposób nie rozwiązuje się problemów. Znaczenie „obejścia” nr 3 z @ michael-e2 w stosunku do wszystkich innych jego obejść można odkryć na nowo, kończąc ścieżkę do pełnego rozwiązania elips i przyjmując jako kompletne rozwiązanie i matematykę dotyczącą prawdziwego rozwiązania i połowy - rozwiązanie, które wszystkie inne oferują tutaj. Robienie tego ręcznie to ciężka praca i dużo pisania. W ten sposób proces Mathematica nie kończy zadania matematycznego kompletnym i poprawnym. Po prostu nie śledzi pracy Reduce.

Ale zachowaj jako kwintesencję odpowiedzi, unikaj korzeni w wynikach z Mathematica w większości przypadków w taki sposób, aby nie pojawiały się w Twojej odpowiedzi, jest blisko poprawnego rozwiązania. Dlatego rozsądne może być leczenie wReduce $y$ i $y'$jako niezależne i wprowadź je odpowiednio. Nie ma wbudowanej funkcji przenoszenia pracy wykonanej Reduceza Ciebie na rozwiązaniu z wyjścia Mathematica. Jest to kwestia doświadczenia, które może osiągnąć każdy matematyk. Jak pokazuje odpowiedź @ michael-e2 , może to prowadzić do nowych gałęzi rozwiązań mieszających wszystkie oznaki korzeni. Dlatego ostateczne rozwiązanie jest jedyne w swoim rodzaju, ponieważ nie ma ambiwalentnego znaku przed korzeniami.