%% 24   Lutningsfält
%% 24.1   Introduktion
% En första ordningens differentialekvation är en ekvation som involverar en obekant funktion  y(t)  och dess derivata  y'(t). Att lösa ekvationen betyder att bestämma alla funktioner som  y(t)  kan vara för att ekvationen ska stämma. Normalt går det inte lösa såna ekvationer analytiskt.
% Vi ska här titta på specialfallet när ekvationen är av formen
% y'(t) = F(t ,y(t))\,,
% där  F( t, y)  är ett känt uttryck som inolverar  t  och  y\,. Ekvationen betyder att funktionsvärden till  y'  kan beräknas direkt från funktionsvärden till y. Om t.ex.  y(2) = 3  så är  y'(2) = F(2,3)  om  y(t)  löser differentialekvationen.
% I detta fall kan man ge en geometrisk tolkning av vad en lösning betyder och dessutom ett sätt, som kallas Eulers metod,  att approximativt bestämma funktionsvärden för lösningar.
%
%% 24.2   Lutningsfält till y'=F(t,y)
% Takaways
% Lutningsfält till en differientialekvation av formen y'(t)=F(t,y(t)).
% Jämviktslösning till en differentialekvation.
% Stabil och instabil jämviktslösning.
% Man kan snickra ihop ett kommando  slope_field  som illusterar lutningsfält.
% Antag att ekvationen är y'= F(t,y)=t^2(y-2)(3-y). Om y  löser den ekvationen och t.ex.  y(1)=4,  så vet också vad y'(1) är, nämligen F(1,4)=1\cdot 1\cdot(-1)=-1. Det betyder att tangenten i  punkten (1,4) på grafen til en lösning ska ha lutning  -1. Uttrycket  F(t,y)  bestämmer alltså vilken lutning grafen till en lösning ska ha i  punkter i  t,\,y-planet .
%
% Ett sätt att tänka på situationen är att  F(t,y)  ger en på förhand bestämd lutning i punkter i  t,\,y-planet  och att hitta en lösning till ekvationen y(t)'=F(t,y(t)) betyder att hitta en funktion y(t), sån att lutningen på grafen till den vid varje tidpunkt  t  stämmer överens med den på förhand givna lutningen F(t,y).
% Genom att välja ut ett rutnät av punkter i t,\, y-planet och i dessa markera ett kort linjestycke med den lutning som F föreskriver kan vi få en geometrisk uppfattning av ekvationen.
% Matlab har (tyvärr) inget färdigt kommando för att illusterar ett lutningsfält, men det går att fixa med ett mer generellt kommando  quiver . För att det inte ska dra iväg förlångt med detta kommer här en färdig modifiering som ger riktningsfält med kommandot  slope_field .
% function slope_field( f , tv , yv )
%     [ Tv , Yv ]=meshgrid( tv , yv ) ; % Skaffar noder
%     t= 1 ./ sqrt( 1 + f( Tv , Yv ).^2 ) ; % Normerar vektorn (1,f(t,y))
%     u= f(Tv,Yv) ./ sqrt(1 + f( Tv , Yv ).^2 ) ; %
%
%     q = quiver( Tv ,Yv , t ,u , 'b') ; % Illusterar halva sträckorna
%     q.AutoScaleFactor = 0.2 ; % Skalar om
%     q.ShowArrowHead = 'off' ; % Inga pilar
%     q.AlignVertexCenters = 'on' ;
%     hold on
%
%     q = quiver(Tv,Yv,-t,-u, 'b') ; % Illusterar andra halvan
%     q.AutoScaleFactor = 0.2 ;
%     q.ShowArrowHead = 'off' ;
%     q.AlignVertexCenters = 'on' ;
%     hold off
% Det är inte meningen att du ska förstå koden, men spara den som  filen  slope_field.m  och se till att den ligger i Matlabs sökväg. Kommandot
% slpoe_field( F , tv , yv ) ,
% där  F  är en anonym funktion som beror på  t  och  y ,  tv  och  yv  radvektorer med t- respektive y-värden, ger nu en illustration av lutningsfältet till differentialekvationen y'(t) = F(t,\, y(t)).
%
% Uppgift 24. 1
% a) Illustera lutningsfältet till differentialekvationen  y'=t^2(y-2)(3-y)
% över området 0\leq t\leq 3,\,\,-1\leq y\leq 4.
% b) Lösningarna till ekvationen kan visas vara
% y(t) = \frac{3C\mbox{e}^{t^3/3}+2}{C\mbox{e}^{t^3/3}+1},
% där  C  är en godtycklig konstant, samt lösningen att y  är konstant  3. Plotta i samma figur de lösningar där y(0)=3/2,\, y(0)= 5/2, samt y(0)= 7/2.
% Lösning



% Lösningar till en differentialekvation som är konstanta kallas jämviktslösningar. De är oföränderliga över tiden (om vi nu använder den som variabel). I uppgift 23.1 fanns det två jämviktslöningar: y(t)=2 samt y(t)=3. Lösningen y(t)=3 har egenskapen att varje lösning som nån gång ligger i närheten av den kommer att ha en y(t)=3 som asymptot i \infty. En sån jämviktslösning kallas stabil. Jämviktslösningen y(t)=2 är instabil.
%
%% 24.3   Eulers metod
% Takeaways
% Metoden ger approximativa funktionsvärden till lösningar till y'=F(t,y).
% Man startar i en punkt (t_0,y_0) som man väljer. Den lösning man då vill approximera är den där y(t_0)=y_0. Punkten (t_0,y_0) kallas lösningens begynnelsevärde.
% Man bestämmer en steglängd \Delta t och flyttar sig successivt i t-led med detta tillskott, samtidigt som man använder riktningsfältet F(t,y) för att flytta sig i y-led.
% Kortare steglängd ger bättre approximtion av en analytisk lösning.
% Låt oss säga att vi vill lösa differentialekvationen  y'(t)=F(t,y(t)),  och söker en lösning sån att vid en fix tidpunkt  t_0  är  y(t_0)=y_0  (där  y_0 är på förhand givet). Vi kan nöja oss med approximativa funktionsvärden till en sån lösning när t\geq t_0.
% Vi startar alltså i punkten (t_0,y_0) på grafen till vår lösing. Där har vi den på förhand givna lutningen F(t_0,y_0), som lösningen ska ha som derivata. Vi bestämmer oss för att följa linjen genom (t_0,y_0) under en kort tid \Delta t (steglängd) och kommer då till punkten (t_0+\Delta t , y_0+\Delta t\,F(t_0,y_0)) som vi skriver som (t_1, y_1). Eftersom vi följer en tänkt tangentlinje till lösningen under en kort tid bör  y_1 inte vara så längt från y(t_1). Denna idé brukar kallas Eulers metod och kan upprepas, så att vi får (t_2, y_2) från (t_1,y_1) och så vidare.
%
% Vi har att
% \begin{array}{rcll}
% t_0 & y_0 & \mbox{ är givna } & \\
% t_i &=& t_{i-1}+\Delta t & \mbox{ när }\, i\geq 1\\
% y_i & =& y_{i-1} + \Delta t \, F(t_{i-1}, y_{i-1}) & \mbox{ när }\, i\geq 1\\
% \end{array}\
% Om \Delta t är litet bör y_i\approx y(t_i) och vi kan illustrera grafen till y(t) utan att lösa differentialekvationen.
%
% Uppgift 24. 2
% Bestäm, med Eulers metod, approximativt de lösningar till  y'(t)=t(y+2)(y-3)/10,  som uppfyller att
% a) y(0)= 2.7,  0\leq t \leq 6,
% b)  y(-1) = 3.1, 0\leq t \leq 5.
% Rita graferna! Använd två olika steglängder i varje deluppgift.
% Lösning
%
% OBS! I videon görs ett misstag i for-satsen för Eulers metod. Det ska (t.ex) vara
% for t = tv1( 1:end-1 )
% och inte  for t = tv1( 2:end-1 ). Misstaget upprepas flera gånger.



%
%% 24.4    ODE45
% Takeaways
% Kommandot  ode45( F , [a , b ], y_0) levererar en radvektor med tidpunkter mellan  a och b och en radvektor med approximativa funktionsvärden i dessa tidpunkter för den lösning till y'=F(t,y) som har y(a)=y_0 .
% F  ska vara en anonym funktion som beror på  t  och  y (eller andra namn på de två variablerna).
% Matlab har förstås (flera) inbygda approximativa lösare av differentialekvationer, eftersom de är så viktiga i tillämpningar av matematik.
% En av dem ges av kommandot
% [ tv , yv ] = ode45( F , [a, b], y_0)
% Här ska  F  vara en anonym funktion som beror av t och y ,  [a, b]  en radvektor som specar start- och sluttid (dvs tidsintervallet) och  y_0  lösningens värde tid tiden t=a.
%
% Uppgift 24.4    Lös uppgift 24.3 a)  igen och lägg till grafen för den lösning som  ode45 ger.
% Lösning
% OBS! I videon görs ett misstag i for-satsen för Eulers metod. Det ska (t.ex) vara
% for t = tv1( 1:end-1 )
% och inte  for t = tv1( 2:end-1 ). Misstaget upprepas flera gånger.



%
%
%
%
%