2013-06-14 10 views
7

Próbuję użyć ode45 rozwiązać system ODE to:Matlab: Czy możliwe jest numeryczne rozwiązanie układu odów z mieszaniną warunków początkowych i końcowych?

[X,Y]= ode45(@sys,[0, T],y0); 

gdzie

function dy = sys(t,y) 

     dy(1) = f_1(y) 
     dy(2) = f_2(y) 
     dy(3) = f_3(y) 
end 

Problemem jest to, że funkcja ode45 wymaga y0 być wartościami początkowymi [y_1(0), y_2(0), y_3(0)], natomiast w moim system, mam dostępne tylko wartości [y_2(0), y_3(0), y_3(T)].

Matematycznie, ten zestaw warunków początkowych/końcowych powinien wystarczyć do ustalenia systemu, ale czy jest jakiś sposób, w jaki mogę z nim pracować przez ode45 lub jakąkolwiek inną funkcję w MATLAB?

Dzięki!

+0

Naprawdę jestem zainteresowany tym pytaniem, ale obawiam się, że nie mogę pomóc; Nigdy nie spotkałem się z tego rodzaju problemem przed ... Wiem, że 'ode45' może integrować się wstecz (po prostu użyj' tspan = [tend tstart] '), więc możesz wymyślić iteracyjny schemat, aby uzyskać' y_1' w taki sposób, że 'y_3 (0) 'i' y_3 (T) 'są spełnione. Nie trzeba dodawać, że będzie to * bardzo * powolne i dość niezgrabne, ale będzie to ** rozwiązanie **. Będę miał na to oko :) Czy mógłbyś zamieścić równania i warunki początkowe/końcowe? –

+0

@RodyOldenhuis Myślę, że znalazłem sposób na rozwiązanie tego problemu. – Vokram

Odpowiedz

6

Po kopaniu w dokumentacji Matlab, myślę, że bardziej eleganckim sposobem jest użycie funkcji bvp4c. bvp4c to funkcja zaprojektowana specjalnie do obsługi problemów związanych z wartością graniczną w ten sposób, w przeciwieństwie do ode**, które są tak naprawdę tylko w przypadku problemów z wartością początkową. W rzeczywistości istnieje cały zestaw innych funkcji, takich jak deval i bvpinit w Matlab, które naprawdę ułatwiają korzystanie z bvp4c. Oto link to the Matlab documentation.

wyślę krótki (i może trochę wymyślony) Przykład tutaj:

function [y1, y2, y3] = test(start,T) 

solinit = bvpinit(linspace(0,3,10), [1,1,0]); 
sol = bvp4c(@odefun,@bvpbc,solinit); 

tspan = linspace(start,T,100); 
S = deval(sol, tspan); 
y1 = S(1,:); 
y2 = S(2,:); 
y3 = S(3,:); 

plot (tspan,y1) 

figure 
plot (tspan,y2) 

figure 
plot (tspan,y3) 


%% system definition & BVCs 

function dydx = odefun(t,y) 
    dydx(1) = y(1) + y(2) + t; 

    dydx(2) = 2*y(1) + y(2); 

    dydx(3) = 3 * y(1) - y(2); 
end 

function res = bvpbc(y0,yT) 
    res= [y0(3) yT(2) yT(3)]; 
end 

end 

W test wyjść Funkcja 3 wektory punktów rozwiązania dla y1, y2 i y3 odpowiednio, a także ich działki.

Oto zmienne ścieżki nakreślonego przez Matlab: y1 path y2 path y3 path

Również znalazłem this video by professor Jake Blanchard z WMU być bardzo pomocne.

+0

+2: Wygląda na to, że polizałeś! Nauczyłem się czegoś nowego dzisiaj, dzięki za to :) –

+0

+2 za zdziwienie, że @RodyOldenhuis nie wiedział o solvercie BVP! :-) Możesz zaakceptować własną odpowiedź, jeśli chcesz. – horchler

3

Jednym ze sposobów jest użycie shooting method do iteracyjnego rozwiązywania dla nieznanego stanu początkowego y_1(0) w taki sposób, że pożądany stan końcowy y_3(T) jest spełniony.

Wpływy iteracji rozwiązując równania nieliniowego F = 0:

F(y_1(0)) = Y_3(T) - y_3(T) 

gdzie wielkie funkcja Y jest roztwór otrzymany przez integrację Ody na przodu w czasie od zbioru warunków początkowych. Zadanie polega na odgadnięciu wartości y_1(0), aby uzyskać F = 0.

Ponieważ obecnie rozwiązujemy równanie nieliniowe, obowiązują wszystkie typowe podejścia. W szczególności, możesz użyć metody bisekcji lub Newtona, aby zaktualizować swoje przypuszczenia dotyczące nieznanego stanu początkowego y_1(0). Uwaga parę rzeczy:

  1. Ody'S zintegrowany na [0,T] w każdej iteracji nieliniowych rozwiązać.
  2. Może być wiele rozwiązań dla F = 0, w zależności od struktury twoich ODE.
  3. Metoda Newtona może zbiegać się szybciej niż bisekcja, ale może również być niestabilna numerycznie, chyba że możesz podać dobry początkowy domysł dla y_1(0).

Przy użyciu istniejących funkcji MATLAB, nielimitowy solver skalarny FMINBND może być dobrym wyborem jako nieliniowy solver.

Mam nadzieję, że to pomoże.

+0

+1: Byłoby to moje podejście, ale myślę, że odpowiedź własna Vokrama powinna być głosowana na górę :) –

Powiązane problemy