1. No category

HH/ITE/BN
Ordinära differentialekvationer och Mathematica
Något om (ODE) och Mathematica
Bertil Nilsson
2016-01-01
1
2
Ordinära differentialekvationer och Mathematica
HH/ITE/BN
ť Förord
På följande sidor presenteras en elementär "streetwise guide" till ordinära differentialekvationer med flitig användning av Mathematica. Framställningen är fåordig, fri från pedanteri men i någon mening fullständig. Det man väsentligen behöver veta om begrepp,
terminologi, beteckningar och teori för att modellera och lösa problem i framtida kurser och yrkesliv som ingenjör, naturvetare eller
lärare klarläggs och typiska exempel ges.
ť Introduktion och terminologi
En differentialekvation (DE) är en ekvation som innehåller en okänd (sökt) funktion och dess derivator. Uppgiften är att bestämma
funktionen så att (DE) uppfylls. Naturligtvis kan man också ha system av differentialekvationer, där man söker flera obekanta
funktioner likt variabler i ett vanligt ekvationssystem.
Många frågeställningar inom naturvetenskap och teknik leder naturligt till modeller som innehåller differentialekvationer. Utan
dessa skulle man behöva genomföra mängder med laboratorieförsök och all slags planering och förutsägelser bli tids- och arbetskrävande. För att kunna planera behöver vi kunna förutsäga saker och ting. I vissa fall är det omöjligt att vänta. Exempel på en fråga
som inte kan vänta till imorgon är; Hur blir vädret imorgon? Att använda (DE) som redskap för att härma verkliga förlopp brukar
kallas simulering och därmed tillhör de ett av de allra viktigaste matematiska hjälpmedlen för en ingenjör och för produktutvecklande företag handlar det om avgörande konkurrenskraft. Därför är det viktigt att lära sig
ställa upp en differentialekvation som modell för en fysikalisk situation,
finna analytiska lösningar till några vanliga enkla typer av differentialekvationer,
använda Mathematica för att bestämma analytiska lösningar och numeriska lösningar då differentialekvationen är för svår eller
omöjlig att lösa analytiskt. Det sistnämnda är naturligtvis det som är vanligast i verkliga livet,
tolka och utvärdera lösningarna i förhållande till den ursprungliga fysikaliska situationen.
Om y x är den sökta funktionen i en (DE) kallas x för den oberoende variabeln. Ofta används tiden t som oberoende variabel
eftersom många (DE) gestaltar hastighetsförlopp och då är ju derivator med avseende på tiden inblandade. Beror y av endast en
variabel säges (DE) vara en ordinär differentialekvation (ODE). I problemformuleringar brukar man friskt blanda de olika
skrivsätten för derivator. Likaså brukar det inte heller vålla någon förvirring att utelämna den oberoende variabeln.
Exempel: Samma (ODE) i olika skepnader
y' x
yx
x,
y'
y
x,
y
x
y
x
Vi söker funktionen y x där x är den oberoende variabeln.
Om y beror av flera variabler säges (DE) vara en partiell differentialekvation (PDE). Dessa är mycket vanliga i ingenjörstillämpningar, t.ex. hållfasthetsberäkningar med den så kallade Finita elementmetoden. Som namnet (PDE) antyder har vi att göra med
partiella derivator. De faller helt utanför ramen för denna kurs.
Exempel: Vi söker utböjningen u hos en vibrerande sträng på en gitarr eller harpa.
Uppenbarligen varierar den både med läget x längs strängen och tiden t, alltså u x, t .
Den klassiska PDE som beskriver strängens rörelse i rummet och tiden kallas
för vågekvationen
2u
t2
c2
2u
x2
där c är en positiv konstant som bär karakteristiska data för strängen,
så som längd, diameter och elasticitet.
En (DE) säges vara av ordning n om n:te derivatan är den högsta förekommande derivatan av y. Man talar om första, andra osv.
ordningens differentialekvation.
Exempel: Vi har exempelvis att y' x
y'' x
x
yx
yx
yn x
är en första ordningens ODE
är en andra ordningens ODE
x
y '' x
sin x
är en n:te ordningens ODE
En (ODE) säges vara linjär om den kan skrivas på formen
n
i 0 ci
x yi x
f x
c0 x y x
c1 x y' x
där ci x och f x är funktioner av x. Annars kallas den olinjär. Då f x
cn x y n x
f x
0 säger vi att (ODE) är homogen annars inhomogen.
HH/ITE/BN
Ordinära differentialekvationer och Mathematica
Exempel: Tydligen är 
x2 y '' x
c2 x
sin x y ' x
c1 x
3
tan x en linjär (ODE) av andra ordningen. Även de tre differentialekvationerna i
f x
föregående exempel är linjära.
Exempel: Vi ser att y ' x y2 x
y '' x y x y x
y '' x sin y x
x
är en första ordningens olinjär ODE
x är en andra ordningens olinjär ODE
0 är en homogen andra ordningens olinjär ODE
Derivator med avseende på tiden är mycket vanliga i fysik. Ofta skrivs de med "prickar", t.ex. Newtons accelerationslag m x F
som är en andra ordningens (ODE). Här är x läget eller koordinaten för kroppen med massan m som utsätts för kraften F i positiv xx
t
riktning. Med definition av hastighet och acceleration har vi prickbeteckningarna x
("x-prick") för hastigheten och x
x
t
2x
t2
("x-prickprick") för accelerationen.
Om en funktion satisfierar en differentialekvation kallas den för en lösning eller partikulärlösning. Exempelvis ser vi efter insättning att y x x3 är en partikulärlösning till y ' x 3 x2 . Men även y x x3 3 och y x x3 7 eller y x x3 C1 , där C1 är en
godtycklig konstant. Mängden av alla partikulärlösningar kallas för den allmänna lösningen eller lösningsskaran. Att lösa en
differentialekvation är att finna den allmänna lösningen. I detta fall är alltså y x x3 C1 den allmänna lösningen till den linjära
ordinära differentialekvationen y' x 3 x2 och representerar genom olika val av C1 ett oändligt antal partikulärlösningar. I figuren
nedan ses några av dem.
yx
30
I den allmänna lösningen till en differentialekvation
ingår alltid lika många godtyckliga konstanter Ci som
differentialekvationen har ordning.
3
2
1
20 x 3 20
x3 15
10 x 3 10
x3 5
x3
x3 5 1
10 x3 10
x 3 15
20 x3 20
2
3
x
30
Ibland har vi också en så kallad trivial lösning y x 0. Exempelvis har differentialekvationen y ' x y x 0 den triviala och
(oftast) för ingenjörer ointressanta lösningen y x 0, vid sidan om den mer intressanta y x C1 x . Att båda är lösninger inses
efter insättning.
ť Begynnelsevärde och randvillkor
För att få den speciella lösning som motsvarar lösningen till ett fysikaliskt problem måste vi fixera samtliga konstanter Ci genom att
vid någon känd ögonblicksbild applicera så kallade begynnelsevärden (BV) eller begynnelsevillkor. Namnet kommer sig av att
man oftast känner tillståndet då man börjar studera systemet, det vill säga i begynnelsen. För entydighet måste vi ange lika många
begynnelsevärden som vi har ordning. En uppsättning ODE
BV kallas för ett begynnelsevärdesproblem (BVP). Ibland fixeras
konstanterna vid olika lägen eller tidpunkter, då talar man om randvärden (RV) och randvärdesproblem (RVP). Ofta ser man en
kombination av dem och ibland innehåller (ODE) ytterligare konstanter som ska fixeras med hjälp av kända tillstånd under resans
gång. I brist på fantasi kallas även dessa tillstånd lite vanvördigt för (RV).
Exempel: En bil påverkas av en konstant framdrivande kraft Fd ,
Fa
luftmotståndet som är proportionellt mot farten i kvadrat Fa cx2
och rullmotståndet som är proportionellt mot farten Fr kx. Sök
det begynnelsevärdesproblem BVP som bestämmer bilens läge
som funktion av tiden, det vill säga x t .
m
Fr
Lösningsförslag: Skådespelet vi söker beskrivs av Newtons accelerationslag mx
x
t
x hastigheten,
2x
t2
Fd
x
F, där x t är läget i ett jordfast koordinatsystem,
x accelerationen och F summan av de krafter som verkar i koordinatriktningen x. Eftersom det är en andra
ordningens differentialekvation så vi behöver två (BV) för entydig lösning. Antag därför att bilen startar från vila och att vi rullar ut
måttbandet samtidigt som vi startar klockan. Då får vi
BVP
mx
x0
Fd
0
x0
0
cx2
kx
ODE
BV
.
4
Ordinära differentialekvationer och Mathematica
HH/ITE/BN
Vi återkommer med mer konkreta exempel på hur Ci bestäms med hjälp av (BV) och (RV) under varje avsnitt längre fram.
Notera att begynnelsevärden (BV) och randvillkor (RV) skall appliceras på den allmänna lösningen till en
differentialekvation. För att få en entydig lösning behövs det lika många BV
RV som vi har Ci ,
det vill säga lika många som differentialekvationen har ordning. Man börjar alltså med att bestämma
den allmänna lösningen till (ODE) därefter appliceras BV
RV .
ť Ordinära differentialekvationer och Mathematica
En stor klass av ingenjörsproblem kan modelleras av så kallade separabla första ordningens (ODE), linjära första ordningens (ODE)
eller linjära andra ordningens (ODE) med konstanta koefficienter. Men innan vi ger oss i kast med dessa och en uppsjö exempel kan
det vara läge att se vad en av våra bästa vänner Mathematica har att säga i ärendet.
Vi börjar med en enkel första ordningens differentialekvation y' x
DSolve y ' x
yx
y x
yx
0 och dess lösning.
0, y x , x
x
c1
I Mathematica används funktionen DSolve för att lösa en stor klass av differentialekvationer, allt från enkla separabla och linjära
av godtycklig ordning till mycket komplicerade olinjära. Den löser även system av differentialekvationer och är mycket lättanvänd.
Strängt taget handlar det om att skriva av rätt! Lägg märke till att Mathematica förstår den vanliga nomenklaturen med ' när vi menar
derivata. Notera de nödvändiga []-parenteserna eftersom y x ska vara en funktion av x! För övrigt kan man naturligtvis använda
vilka namn man vill. Observera dubbla likhetstecken eftersom det är en ekvation! Det är inte bara namnet som antyder släktskap
med Solve, utan även hantering av indata och resultat. Som vanligt gäller att när man väl förstått filosofin bakom Mathematica är
det mesta självklart! Här kommer några till. Notera speciellt, som sig bör, att vi får lika många obestämda konstanter ci , vilka vi
kallar Ci vid handräkning, som vi har ordning. Det är ordning på Mathematica
DSolve y ' x
y x
sin x
cos x 
2
x2  y ' x
DSolve1
y x
Sin x , y x , x
1
x
c1
y x
1
c1 x2
1
c2
2t
1 tan
1
1
x2  ArcTan x , y x , x
x 2 
2
DSolvex '' t
x t
x2
2xy x
t
4 x' t
2t
c1
1
2t
4x t
, x t , t
2t 2
t 
2
Man kan givetvis ta med begynnelsevärden eller randvärden för att få konstanterna bestämda. Dessa paketeras då tillsammans med
differentialekvationen i en lista. Tänk på att även begynnelsevärdena ska anges som ekvationer, det vill säga med två likhetstecken.
Så begynnelsevärdesproblemet
y' x y x sin x ODE
BVP
y0 2
BV
är bara att skicka rakt in i Mathematica och piggar sedan upp oss med en bild över situationen.
yAvx DSolve y ' x
y x
Sin x , y 0
2 ,y x ,x
Plot y x
. yAvx, x, 0, 5 , PlotStyle Red, AxesLabel
1
y x
5
x
sin x
cos x 
2
yA x
2.0
1.5
1.0
0.5
1
0.5
2
3
4
5
x
Simplify
"x", "yA x "
HH/ITE/BN
Ordinära differentialekvationer och Mathematica
5
Då varken vi eller Mathematica kan finna en analytisk lösning finns funktionen NDSolve till vår hjälp för att göra en ren numerisk
lösning. Utdata från denna är en InterpolatingFunction som kan verka lite märkvärdig innan man blivit vän med den. Den
fungerar dock som vilken annan funktion som helst. För övrigt är den ett kraftfullt redskap om man vill göra interpolation i t.ex.
mätdata. Som exempel kör vi en repris på begynnelsevärdesproblemet ovan. Det enda som skiljer i menyn som serveras jämfört med
DSolve är att man, likt Plot, måste ange i vilket intervall man vill att spektaklet ska utspela sig.
NyAvx
NDSolve
y x
y' x
y x
Sin x , y 0
Domain: 0. 5. InterpolatingFunction
Plot y x
2 , y x , x, 0, 5
Output: scalar
. NyAvx, x, 0, 5 , PlotStyle
 x 
Blue, AxesLabel
"x", "yN x "
yN x
2.0
1.5
1.0
0.5
1
2
3
4
5
x
0.5
Den triviala lösningen y x
DSolve y ' x
yx
0 levereras normalt inte av Mathematica. Exempelvis det inledande exemplet i repris
y x
0, y x , x
x
c1
Notera dock att om vi med (BV) startar någonstans på x-axeln, det vill säga y x0 0, så är vi dömda att stanna kvar där enligt den
triviala lösningen. Så återigen, olika (BV) kan ge helt olika lösningskurvor. Detta inser naturligtvis Mathematica.
DSolve
yx
y' x
y x
0, y 17
0 ,y x ,x
0
ť Riktningsfält
En första ordningens (ODE) kan alltid skrivas på den generella implicita formen
g x, y, y'
0
I många fall kan vi lösa ut y' explicit och skriva den på formen
y'
f x, y
Lösningarna till denna (ODE) kommer att utgöras av kurvor i xy-planet, en för varje konstant C1 . Vi vet ingenting om deras utseende
men vi kan dra följande slutsats. Om en lösningskurva går genom en viss punkt x0 , y0 så måste kurvans tangent i den punkten ha
riktningskoefficienten f x0 , y0 . Detta betyder att om vi i ett antal xy-punkter, oftast ett rutnät, markerar lösningskurvans lutning i
dessa punkter får vi ett så kallat riktningsfält. Detta ger ofta en kvalitativt god bild över situationen. Man kan skissa den lösningskurva man får för givna (BV) och man kan studera dramatiken i omgivningen. Ibland gäller nämligen att små störningar i (BV)
ger en helt ny lösningskurva, så kallat kaotiskt beteende. Kanske har du hört talas om "butterfly"-effekt som utgör en av flera
svårigheter med att göra långa väderleksprognoser. Vi lättar upp stämningen med ett
y'
Exempel: Rita riktningsfält och lösningskurva till BVP
y0
xy
1
2
Lösningsförslag: Först lösningen till (BVP)
1
yAvx
DSolvey ' x
xy x , y 0
, y x , x
2
x2
2

y x
2
ODE
BV
.
6
Ordinära differentialekvationer och Mathematica
HH/ITE/BN
Sedan ritar vi in denna tillsammans med riktningsfältet som är färglagt efter riktningskoefficientens storlek.
Show Plot y x
. yAvx, x, 0, 1.05 ,
VectorPlot 1, x y , x, 0, 1 , y, 0, 1 , VectorColorFunction
AspectRatio Automatic, AxesLabel
"x", "y"
Hue ,
y
0.8
0.7
0.6
0.0
0.2
0.4
0.6
0.8
x
1.0
ť Separabel första ordningens (ODE)
En (ODE) kallas separabel om den kan skrivas på formen
y
g y
f x
x
där g y och f x är kända och kontinuerliga funktioner i var sitt intervall. Då har de primitiva funktioner G y respektive F x i
dessa intervall. Med definition av primitiv funktion och kedjeregeln i färskt minne får vi så den allmänna lösningen på följande sätt
y
g y
x
f x
g y
y
f x
x
0
G y
y
y
x
KR:
F x
x
0
x
G y
F x
0
G y
F x
C1
G y
x
Namnet separabel kommer sig av att om vi även betraktar måtten x och y som variabler så kommer x och y att separeras av
likhetstecknet
y
g y
Eftersom G y
g y
y och F x
f x
f x
x
g y
y
f x
x
x får vi i praktiken den allmänna lösningen genom att först separera variablerna och
sedan "hänga på två integraltecken och en konstant i högerledet"
g y
y
f x
x
C1
Samma svar kommer vi fram till genom att formellt integrera båda sidor av g y
y
x
f x med avseende på x och sedan betrakta det
som variabelsubstitution i vänsterledet
g y
y
x
x
f x
x
C1
g y
y
f x
x
C1
Strängt taget får vi en integrationskonstant för varje obestämd integral. Men dessa kan bakas samman till en enda, vilket ju också är
behovet med tanke på att vi har att göra med en första ordningens (ODE). Vi ser också att lösa en separabel (ODE) faller tillbaka på
att kunna integrera. Mot bakrund av hur lösningsmetoden är uppbyggd kommer den allmänna lösningen naturligt ut på implicit form
h x, y 0. Vid handräkning brukar man nöja sig med detta svar huruvida det inte är uppenbart att kunna lösa ut y x explicit, medan
Mathematica däremot tycker om att späka sig till det yttersta för att i varje läge försöka leverera explicit form, inte sällan med hjälp
av diverse exotiska funktioner. Det kan därför ibland vara svårt att jämföra lösningarna. Vi sammanfattar
En separabel differentialekvation g y
x2
Exempel: Lös differentialekvationen y'
4x
Lösningsförslag: Separabel
x2 x
y
4x
x2
0, y x , x
"integreras direkt".
DSolvey ' x
x3
y x
c1
3
4x
2 x2 
y
x
f x har allmänna lösningen g y
y
f x
x
C1
0.
y
2 x2
1
3
x3
C1 . För den här typen av enkla separabla brukar man säga att de
HH/ITE/BN
Ordinära differentialekvationer och Mathematica
Exempel: Lös differentialekvationen y'
12 x2
8 x3
0 med begynnelsevillkoret (BV) y 1
Lösningsförslag: Vi har ett (BVP) och börjar alltid med att lösa (ODE) som är separabel
y
12
1
3
x3
8
1
4
x4
4 x3
C1
2 x4
y x
2 x4
2 x3
12 x2
8 x3
4 13
2: 2
0, y 1
c1
8 x3 x
2 14
C1
C1
4 x3
4. Så lösningen till (BVP) är y
2 x4
4.
2, y x , x
5 x2
cos 3x .
Lösningsförslag: Separabel eller "direkt integration"
y x
12 x2
y
2
Exempel: Lös differentialekvationen y'
DSolvey ' x
2.
C1 . För den här typen av enkla separabla brukar man säga att de "integreras direkt".
Konstanten C1 fixeras sedan av (BV) y 1
DSolvey ' x
7
5 x2
5 x3
1
3
3
y
5x2
cos 3x
y2
1
2
x
y
5
3
x2
y2
x3
1
3
sin 3x
C1 .
Cos 3 x , y x , x
sin 3 x 
x
.
y
Exempel: Lös differentialekvationen y'
Lösningsförslag: Separabel y'
x
y
dem ny innebörd, i detta fall 2C1
C12 är ett mycket vanligt hyss i branschen. Vi har ju att göra med en godtycklig konstant. Om det
y y
1
2
x x
x2
C1
C12 . Att döpa om konstanter eller att ge
är ett (BVP) vi löser kommer det att bli rätt till slut ändå. Mathematica levererar lösningen på explicit form.
x
, y x , x
DSolvey ' x
y x
y x
2 c1
x2 , y x
2 c1
Geometriskt betyder den implicita formen x2
med centrum i origo och valfri radie C1 .
x2 
C12 cirklar
y2
y
x
Exempel: Lös differentialekvationen y'
Lösningsförslag: Separabel y'
x
y
x
.
y
y y
x x
1
2
y2
1
2
x2
en hyperbel där orienteringen bestäms av tecknet på C1 . Fallet C1
nedan. Mathematica levererar lösningen på explicit form.
C1
x2
y2
0 ger strålarna y
C1 . Geometriskt betyder den implicita formen
x som också är hyperbelns asymptoter, se
x
, y x , x
DSolvey ' x
y x
y x
2 c1
x2 , y x
2 c1
x2 
y
y
y
x
C1
0
x
C1
0
x
C1
0
8
Ordinära differentialekvationer och Mathematica
2y
Exempel: Lös differentialekvationen x3
Lösningsförslag: Separabel x3
x
y
2
2 C1
2 C1
y
x
y'
.
2y
y'
1
y'
2y
x3
2

2 C1
HH/ITE/BN
1
2y
C1 
y
1
y
x
C1
2
1
2
x
x3
1
ln y
2 x2
C1
ln y
x2
C1
1
2 C1
x2
.
2y x
DSolvex3
, y x , x
y' x
1
y x
c1
x2

Exempel: Lös differentialekvationen y'
Lösningsförslag: Separabel y'
DSolve y ' x
2x
2x
2x
xy.
y'
xy
1
x
y 2
y
y 2
x x
ln y
1
2
2
x2
y
C1
2.
2
xy x , y x , x
x2
y x
c1
2
2
Exempel: Lös differentialekvationen yy'
Lösningsförslag: Separabel yy '
x
x
y'
y'.
y
1 y'
x
y
1 y
1
2
x x
y2
1
2
y
x2
y2
C1
2y
x2
C1 .
Implicit form räcker gott! Men med kvadratkomplettering av y får vi y2 2 y 12 12
x2 C1
x2
y 1 2 C1 1
2
2
2
2
2
2
x
y 1
C1 vilket vi identifierar med cirkelns ekvation x x0
y y0
r , det vill säga lösningskurvorna är cirklar med
centrum i 0, 1 , se nedan. Mathematica levererar lösningen på explicit form.
DSolve y x y ' x
y x
x2
1
x
2 c1
y' x , y x , x
x2
1 , y x
2 c1
1
1
y
2.0
1.5
1.0
0.5
1.0
x
Exempel: Lös differentialekvationen y'
x
Lösningsförslag: Separabel, ty vi har y'
y
1
y
x
y
x
u
1
y
x y
osv.
0.5
1.0
x
.
x y
ln 1
0.5
x y
y'
y
x
x
C1
x
y'
ln 1
y
y
x
1
1
y
1
x
y
C1 . Implicit form räcker!
Mathematica levererar lösningen på explicit form.
x
DSolvey ' x
y x
c1
log1
x
x y x
, y x , x
Simplify

Exempel: Lös differentialekvationen x y ' cos y
sin y
0.
Lösningsförslag: Separabel, ty vi har att x y ' cos y
sin y
C1
ln C1
u
sin y osv.
ln sin y
ln x
C1
0
x y ' cos y
sin y
Mathematica levererar lösningen på explicit form.
DSolve x y ' x Cos y x
y x
sin
1
c1
x 
Sin y x
0, y x , x
sin y
x
cos y
sin y
y
C1 x. Implicit form räcker!
1
x
x
HH/ITE/BN
Ordinära differentialekvationer och Mathematica
Exempel: Lös differentialekvationen x2
1 y '
Lösningsförslag: Separabel, eftersom x2
2 C1
ln C1
ln
y2
x2
1

x2
x y.
1 y '
y2
ln C1
1
y
xy
y2
C1
1
9
x
y
C1 x2
x
x2 1
x2
u
1 osv.
1
lnx2
2
ln y
1 C1
1, Geometriskt betyder detta en kurvskara i xy-planet där varje
bestämt värde på C1 ger en medlem av denna kurvskara. För C1
0 fås hyperbler, C1
0 ellipser och för C1
0 linjen y
0.
Mathematica levererar lösningen på explicit form.
DSolvex2
y x
1 y ' x
x2
c1
x y x , y x , x
1 
x
y' 4
y0 3
Exempel: Lös BVP
ODE
.
BV
Lösningsförslag: Separabel, ty x y' 4
y' 4
3
4 0 C1 C1 7. Så lösningen till (BVP) y
x
DSolve
yx
7
y' x
4, y 0
x
y
4
x
4
7.
x
3 ,y x ,x
C1 . Till slut fixeras så konstanten C1 av (BV) y 0
3:
Simplify
x
4
För att applicera (BV) har vi i detta och tidigare exempel använt oss av obestämd integral. Man bör även bemästra att genomföra
kalkylen med bestämd integral. Egentligen är det samma sak man gör men de olika betraktelsesätten är en fråga om smak och de
dyker upp i diverse olika framställningar. Vi tar dem parallellt
BV med obestämd integral
x
y'
y
x
4
y
y0
4
x
4
3:3
C1
y
4
7
4
x3  y '
1
Exempel: Lös BVP
x
y1
2
BV med bestämd integral
x
Separera
x
y
3
Integrera
y'
4
y
x
4
0
y
3
4
Separera
x
x
Integrera med BV direkt
Gränserna är ögonblicksbilder
Känd respektive allmän
C1
0
C1
BV
y
7
y
Lösningen
till BVP
x2 y
ODE
BV
3
y
x x
0
4
3
y
Insättningstecken, döljer C1 .
x
0
x
4
4
Nu är det bara att hyfsa
4
7
x
Lösningen till BVP
.
Lösningsförslag: Innan vi exercerar våra nyvunna kunskaper ännu en gång, kontrollerar vi facit med lösningen på explicit form.
x3  y ' x
DSolve1
3
22 3
y x
x3
x2 y x , y 1
2, y x , x
1 
BV med obestämd integral
x3  y '
1
1
y
x
y
1
3
ln y
3 ln y
x2 y
ln 1
x
x3
C1
3
3 C1
x
lnab 
y
y1
C1 1
3
2:2
C1
3
y
ln b
3
x
C1 1
4
4 1
ln y
b ln a ,
3
x
1
3
BV
ln y
1
3
Separera
2
x x
1 1 x3
y
y
2
1
3
x
x3
ln 1
Integrera med BV direkt
Gränserna är ögonblicksbilder
Känd respektive allmän
x
1
Insättningstecken, döljer C1 .
ln 2
x3
ln 1
3 ln y
ln 1
ln y3
Lösningen till BVP
x2 y
x  y'
y1
2 y
Integrera
u 1 x3 osv.
ln ab ln a
3 C1 ln C1
3
1
Separera
2
1 x3
ln 1
BV med bestämd integral
3
Nu är det bara att hyfsa
3 ln 2
x3 ln 2
ln 1
y3
13
ln 1
x3
4 1
ln 22
x3
Lösningen till BVP
10
Ordinära differentialekvationer och Mathematica
HH/ITE/BN
ť Linjär första ordningens (ODE)
Enligt definition ovan kan en första ordningens linjär (ODE) skrivas på formen
y'
px y
qx
där p x och q x är kända och kontinuerliga funktioner. Som lösningsstrategi ska vi välja att multiplicera båda sidor i differentialekvationen med en funktion f x och sedan jämföra vänsterledet med derivatan av produkten f x y. Vi får då
f x y'
f x px y
f x qx
Lika om
f' x
f x px
x
Vi ser att om vi väljer f x så att f ' x
f x y
f x p x blir
x
f' x
f x y
f x q x . Detta är en separabel (ODE) vars allmänna lösning är
 f x qx
x
C1 
f x p x ? Men visst! Detta är också en separabel (ODE) med lösningen
f
f x px
f' x y
1
f x
yx
Frågan är bara om det finns ett f x så att f ' x
f x y'
px
f x
x
ln f x
px
x
C1
px
f x
x C1
Att lösa en linjär första ordningens (ODE) på detta sätt kallas för metoden med integrerande faktor, just på grund av det lämpliga
faktorvalet som gör att lösningen fås efter en integration. Med detta val kallas f för integrerande faktor och vi benämner den med
IF. Vanligtvis väljer man i IF att sätta integrationskonstanten C1 0 eftersom den ändå kan divideras bort, ty
px
IF x
px
x
A y'
x C1
px
px
px y
px
x
y'
x C1
x
px
Aq x
x
A med A
0
Förkorta bort A
px
px y
x
0
qx
Däremot får man aldrig glömma C1 i den allmänna lösningen. Liksom i det separabla fallet gäller det att vara en smidig integrator!
Vi sammanfattar.
En linjär första ordningens differentialekvation y '
Då är y '
px y
qx
IF x y
x
px y
px
q x har integrerande faktorn IF x
IF x q x med allmänna lösningen y x
1
IF x
 IF x q x
x
x
.
C1 .
Lär er inte formeln utan metoden! Glöm inte att multiplicera högerledet q x med IF x före separation!
Exempel: Lös differentialekvationen y'
y
1.
Lösningsförslag: Separabel eller linjär, här väljer vi det senare. Genomför gärna det förstnämnda! Då blir integrerande faktorn
1 x
IF
y 1
C1
x
x
, och därmed
x
x
x
y
1
yx
c1
x
y x
Som vanligt är

x
y
x
C1
1
2x
C1
1
2
y
DSolve y ' x
x
y x
C1
x
x
y
.
1 x
Lösningsförslag: Linjär (ODE) med IF
1
2
x
1, y x , x
Exempel: Lös differentialekvationen y'
y
x
y
.
DSolve y ' x
x
x
Separabel
x
x
, så
x
x
y
x
x
x
Separabel
2x
y
x
.
,y x ,x
x
y x
x
c1

2
Exempel: Lös differentialekvationen y'
5y
x
.
Lösningsförslag: Ej separabel men linjär. Jämför med prototypen y '
x

5x
y
5x
x
Separabel

5x
y
4x
x
5x
px y
y
1
4
qx
4x
y'
C1
x
5 y
y
1
4
x
. Så IF
C1
5x
.
5 x
5x
och
HH/ITE/BN
Ordinära differentialekvationer och Mathematica
DSolve y ' x
x
5y x
11
,y x ,x
x
y x
5x
c1

4
Exempel: Lös differentialekvationen y'
2y
3.
2 x
Lösningsförslag: Separabel eller linjär, som omväxling väljer vi det senare. Då blir integrerande faktorn IF
x

2x
2x
y
3
Separabel
DSolve y ' x
y x
2y x

y
2x
3
x

x
(BV) y 2
x
DSolve
x y'
y2
2y
1
lnx2 
2
3x
1
C1
22
x y' x
x2
x2 , så att
C1
1
2
2
y x
y
y
C1
, så
.

2 ln x
1: 1
2
x
C1
2x
3, y x , x
ODE
.
BV
Lösningsförslag: Efter division med x får vi en linjär (ODE), y '
IF
y
2x
2x
2
Exempel: Lös BVP
2
2x
3
2
3
2x
c1
2x
3
2
2y x
x
x2 y
1

x
x2 3
1
,
x
3
och därmed är vi i välkänd terräng. Vi får här
x2 y
Separabel
x3
1
2
x2
C1
y
1
2
x
C1
x2
. Slutligen med
2.
3x
1, y 2
1 ,y x ,x
Simplify
1
x

2
ť Linjär andra ordningens (ODE) med konstanta koefficienter
Vi ska studera en i ingenjörssammanhang vanligt förekommande variant av andra ordningens linjära (ODE), nämligen då koefficienterna framför y'', y ' och y är reella konstanter. Då så inte är fallet har vi i allmänhet stora bekymmer att vänta och överlåter verksamheten till Mathematica. Alltså
y''
a y'
by
f x
där a och b är reella konstanter. Den allmänna lösningen till denna (ODE) består av två delar
yh x
yx
yp x
där homogena lösningen yh x är lösningen till motsvarande homogena differentialekvation y''
partikulärlösning till den orörda y '' a y' b y
göra en korrekt ansats av partikulärlösningen.
ay'
by
0 och y p x är en
f x . Man börjar alltid med den homogena lösningen eftersom den behövs för att
Homogena lösningen yh x
Studera den homogena differentialekvationen
y ''
ay'
by
0
och byt y mot r och derivator mot potenser, så får vi en vanlig andragradsekvation i r som kallas för den karakteristiska ekvationen. Denna har alltid två rötter som bestäms med välkänd formel
y ''
ay'
by
ar
r
r1,2
r2
0
2
a
2
ar1
b
br0
0
a 2
2
b
Den homogena lösningen beror nu på de två rötterna r1 och r2 . Vi sammanfattar
0
12
Ordinära differentialekvationer och Mathematica
HH/ITE/BN
Till den homogena differentialekvationen
y '' a y ' b y 0
hör karakteristiska ekvationen r2
ar
b
a 2
a
2
0 med rötterna r1,2
2
b.
Den homogena lösningen yh x får vi då som ett av de tre fallen
Fall 1: r1 r2 reella yh x C1 r1 x C2 r2 x
Fall 2: r1 r2 reella yh x
C1 x C2 r1 x
Fall 3: r1,2 Α
Β, Β 0 komplexkonjugerade
Αx
yh x
C1 cos Βx
C2 sin Βx
Lär dessa utantill!
Vi ser att det är den homogena lösningen som bär de två konstanterna C1 och C2 . För att övertyga sig om riktigheten i de tre fallen
och att det inte finns några andra kan vi argumentera på följande sätt.
Gör variabelsubstitutionen y x
Insättning i y''
a y'
by
r2 x
zx
. Derivera produkten en och två gånger. Efter hyfsning har vi
y' x
r2 x
z' x
y '' x
r2 x
z'' x
r2 z x
r22 z x 
2 r2 z ' x
0 ger
r2 x
z '' x
2 r2
r22
a z' x
ar2
b z x 
0
Eftersom r2 x 0 för alla x och r2 är rot till karakteristiska ekvationen kan vi tillsammans med sambandet r1
och koefficienter i en andragradsekvation förenkla
r2 x z'' x

2 r2
0
r22
a z' x
ar2
r2 r1
b z x 
r2
a mellan rötter
0
0
till en linjär (ODE) av första ordningen i z' x .
z' x '
r2
r1 z' x
0
Denna har då enligt tidigare avsnitt lösningen
z' x
som är separabel. Om nu r1
zx
r2 x
r
C1
r2
r1 r2 x
1
Om däremot r1
r2 blir z' x
C1
r1 r2
C2 
r2 x
r1 x
C1 x
r1 r2 x
r
C1
r2
1
C2
C1 
C1
r1 x
C2
r2 x
vilket är fall 1.
C1 och efter separation
zx
varav y x z x
identitet ger
r1 r2 x
r2 får vi
zx
Substituera tillbaka y x
C1
C2
yx
r1 x
C1 x
C2
som är fall 2. Fall 3 är något mer komplicerat men utveckling av fall 1 med hjälp av Eulers
C1
Αx
Αx
Αx
r1 x
C2 r2 x C1 Α Β x C2 Α Β x
C1 cos Βx
sin Βx
C2 cos Βx
C1 C2 cos Βx
C1 C2 sin Βx
Acos Βx Bsin Βx
Eulers identitet
sin Βx
där vi infört nya komplexa konstanter A C1 C2 och B
C1 C2 . Eftersom reella lösningar är de enda fysikaliskt relevanta bör
dessa väljas rent reella och fall 3 skall användas då vi har komplexkonjugerade rötter till karakteristiska ekvationen. De tre möjliga
fallen av homogena lösningar illustreras kvalitativt i följande figurer. Vokabuläret är hämtat från mekanikens studium av svängningsproblem
HH/ITE/BN
Ordinära differentialekvationer och Mathematica
yx
yx
13
yx
x
x
x
Fall 2: Kritisk dämpning.
Lösningskurvan skär x axeln högst
en gång. Typiskt utseende för en ny
hjulupphängning på en bil.
Fall 1: Överkritisk dämpning.
Inte så vanligt förekommande i
ingenjörssammanhang.
Exempel: Sök homogena lösningen till y''
4 y'
5y
Fall 3: Underkritisk dämpning.
Mycket vanligt insvängningsförlopp.
Typiskt utseende för en gammal
hjulupphängning på en bil.
0.
Lösningsförslag: Vi får karakteristiska ekvationen och dess två rötter
r2
Solver2
r
4r
1, r
5
4r
5
0
5
2
r1
r2
3
5
1
5
DSolve y '' x
x
c1
4 2

2

0
Två olika reella rötter, så Fall 1 ger direkt yh x
y x
4
2
r1,2
c2
4 y' x
5x
5y x
5x
C1
x
C2
. Detta upptäcker naturligtvis även Mathematica.
0, y x , x

Exempel: Sök homogena lösningen till y''
4 y'
4y
0.
Lösningsförslag: Vi får karakteristiska ekvationen och dess två rötter
r2
Solver2
r
2, r
4r
4
4r
4
2
4
2
0
r1
r2
2
2
2
DSolve y '' x
2x
4 2
4
2
r1,2
0
Reell dubbelrot, så Fall 2 ger direkt yh x
y x
0
c2 x
4 y' x
C1 x
4y x
C2
2x
. Detta upptäcker naturligtvis även Mathematica.
0, y x , x
Simplify
c1 
Exempel: Sök homogena lösningen till y''
4 y'
13 y
0.
Lösningsförslag: Vi får karakteristiska ekvationen och dess två rötter
r2
Solver2
r
2
4r
13
3 , r
2
4r
13
0
r1,2
4
2
y x
2x
2
3
r1
r2
2
2
3
3
3
4 y' x
c1 sin 3 x
13
0
Komplexkonjugerade rötter, så Fall 3 ger direkt yh x
Mathematica.
DSolve y '' x
4 2
2
13 y x
c2 cos 3 x 
0, y x , x
2x
C1 cos 3x
Simplify
C2 sin 3x . Detta upptäcker naturligtvis även
14
Ordinära differentialekvationer och Mathematica
HH/ITE/BN
Partikulärlösningen y p x
Studera nu den ursprungliga inhomogena differentialekvationen
y''
a y'
by
f x
där den homogena lösningen yh x är bestämd enligt föregående avsnitt. Vi ska nu bestämma en partikulärlösning y p x som beror
av högerledet. Det räcker att hitta en och för enkelhets skull nöjer vi oss här med att behandla några enkla och i praktiken vanligt
förekommande funktioner f x . Övriga (i praktiken alla) fall överlåter vi till Mathematica. Vi sammanfattar resan i en kokbok och
illustrerar sedan med några exempel.
Steg 0: Innan man börjar så kan man dra lite nytta av lineariteten hos (ODE) för att bespara sig alltför volymiösa handräkningar.
Nämligen om f x
f1 x
f2 x
där fi x är funktioner som vi befattar oss med enligt Steg 1 nedan, så gäller
f x
f1 x
f2 x
yp x
y p1 x
y p2 x
Man kan alltså genomföra Steg 1-3 nedan för varje fi x och sedan lägga samman y pi x till den slutliga y p x som motsvarar f x .
Steg 1: Vi väljer att behandla högerled f x enligt kolonn ett. Ansätt y p x enligt kolonn två med okända koefficienter A, B,
som
sedan kommer att bestämmas i Steg 3.
f x
yp x
polynom av grad n
Ansätt ett fullt polynom av samma gradtal n,
det vill säga ett som innehåller alla potenser upp till och med n:
A
k
Αx
C x2
Bx
N xn
Αx
Ansätt A
k cos Ωx
Ansätt A cos Ωx
B sin Ωx
k sin Ωx
Ansätt A cos Ωx
B sin Ωx
Ansätt A cos Ωx
B sin Ωx
k1 cos Ωx
k2 sin Ωx
Exempel: Typiska enkla fall som man bör klara för hand. Redan detta är på gränsen
f x
9
f x
5x
f x
f x
2
7x
3
x
2
3
yp x
A
yp x
A
yp x
4x
A
yp x
A
Bx
Bx
Bx
2
Cx
2
Cx
3
f x
5 cos 4x
f x
3
Dx
2x
f x
f x
12 sin 4x
3 cos 4x
8 sin 4x
2x
yp x
A
yp x
A cos 4x
B sin 4x
yp x
A cos 4x
B sin 4x
yp x
A cos 4x
B sin 4x
Steg 2: Kontrollera om någon term i y p x finns med som en term i yh x , i så fall tillför y p x inte någon ny information. Korrigera
genom att multiplicera y p x med x. Om y p x fortfarande finns i yh x så multiplicera med x igen, å igen
Detta problem uppkom-
mer bara när vi har en exponentialfunktion eller trigonometriskt uttryck i högerledet.
Exempel: Ansätt y p x beroende på f x . Undersök sedan om y p x finns med som termer i yh x och korrigera vid behov. Det är ok
om y p x finns med i yh x som en faktor med annan funktion g x .
yh x
C1
5x
C2
2x
& f x
3
2x
yp x
yp x
yh x
C1 x
C2
2x
& f x
3
2x
yp x
Ax
2
Ax
yh x
2x
C1 cos 4x
C2 sin 4x
& f x
3
yh x
2x
C1 cos 4x
C2 sin 4x
& f x
3 cos 4x
0x
C1 cos 4x
C2 sin 4x
& f x
2x
A
yp x
2x
3 cos 4x
yh x
2x
Ax
yp x
yh x
2x
A
yh x
yh x
2x
yh x
2x
yh x
yp x
yp x
yp x
yp x
A
2x
yh x
A cos 4x
B sin 4x
yh x
A cos 4x B sin 4x yh x
x A cos 4x B sin 4x
yh x
HH/ITE/BN
Ordinära differentialekvationer och Mathematica
Steg 3: Sätt in ansatsen y p x i (ODE) och bestäm de okända koefficienterna A, B,
15
genom att identifiera koefficienter. Detta
innebär att bestämma A, B,
så att vi har lika många "äpplen och päron" på båda sidor om likhetstecknet. Det vill säga för att en
likhet skall gälla för alla värden på den oberoende variabeln måste vi ha exakt samma antal av varje förekommande funktion på båda
sidor om likhetstecknet.
Exempel: Sök y p x då y ''
6 y'
13 y
5x
2 är given.
Lösningsförslag: Först bestämmer vi den homogena lösningen med hjälp av karakteristiska ekvationen
r2
6r
Solver2
r
3
13
0
6 2
6
2
2
13
2 så ansätt enligt tabell y p x
A
6r
13
2 , r
3
r1,2
3
2
3x
yh x
Fall 3
C1 cos 2 x
C2 sin 2 x
0
2
Sedan en partikulärlösning.
Steg 1: Vi har f x
5x
Steg 2: Vi ser att y p x
B x.
yh x .
Steg 3: Identifiera koefficienter. Först bäddar vi upp med derivator
yp x
A
Bx
yp ' x
B
y p '' x
0
Insättning i differentialekvationen
0
6 B
13 A
Bx
5x
2
och identifiering av koefficienter framför varje funktionstyp (FT) i vänsterledet (VL) och högerledet (HL)
FT
VL
HL
1
13 B
5
x
0
6B
x
13 A
2
Detta linjära ekvationssystem ska alltid ge en entydig lösning. Om inte så är ansatsen av y p x felaktig eller felräkning på vägen! Gå
i så fall tillbaka till Steg 1! Här får vi
Solve
13 B
5, 6 B
56
,B

169
56
169
2
5
A
Så y p x
13 A
13
5
13
x varav allmänna lösningen som summan av den homogena lösningen och partikulärlösningen
yh x
yx
3x
yp x
C1 cos 2 x
56
169
C2 sin 2 x
5
13
x
Med Mathematica har vi direkt
DSolve y '' x
y x
c1
3x
6 y' x
sin 2 x
Exempel: Sök y p x då y ''
c2
2 y'
13 y x
3x
5x
2, y x , x
5x
56
13
169
cos 2 x
3y
2 x2
Simplify

5 är given.
Lösningsförslag: Först bestämmer vi den homogena lösningen med hjälp av karakteristiska ekvationen
r2
2r
3
0
r1,2
2
2

2 2

2
3
1
2
r1
r2
3
1
Fall 1
yh x
C1
3x
C2
x
16
Ordinära differentialekvationer och Mathematica
Solver2
r
2r
1, r
3
HH/ITE/BN
0
3
Sedan en partikulärlösning.
2 x2
Steg 1: Vi har f x
Steg 2: Vi ser att y p x
5 så ansätt enligt tabell y p x
A
C x2 .
Bx
yh x .
Steg 3: Identifiera koefficienter. Först bäddar vi upp med derivator
y p x A B x C x2
yp ' x B 2 C x
y p '' x
2C
Insättning i differentialekvationen
2C
2 B
2C x
3 A
Bx
C x2 
2 x2
5
och identifiering av koefficienter framför varje funktionstyp (FT) i vänsterledet (VL) och högerledet (HL)
FT
VL
2
x
HL
3C
1
4C
x
0
x
2C
2
3B
2B
0
3A
5
Detta linjära ekvationssystem ska alltid ge en entydig lösning. Om inte så är ansatsen av y p x felaktig eller felräkning på vägen! Gå
i så fall tillbaka till Steg 1! Här får vi
Solve
3C
17
A
17
27
4C
8
,B
2B
3A
5
3
2
3
x
0, 2 C

9
8
9
3B
2
,C
27
Så y p x
2,
x2 varav allmänna lösningen som summan av den homogena lösningen och partikulärlösningen
yh x
yx
yp x
C1
3x
C2
x
17
27
8
9
x
2
3
x2
Med Mathematica har vi direkt
yAvx
DSolvey '' x
y x
c1
x
3x
c2
2 y' x
2 x2
3y x
2 x2
8x
17
3
9
27
5, y x , x
Simplify

Mathematica har ingen funktion som levererar partikulärlösningen ensam. Detta är av visst intresse ibland, se vidare exempel i
avsnittet "Blandade exempel för resten". Enklaste sättet är i så fall att bestämma allmänna lösningen till (ODE) och sedan släcka ner
den homogena lösningen genom att sätta c1 och c2 till noll.
yAvx . C 1
0, C 2
2 x2
8x
17
3
9
27
y x
0

Även då man räknar för hand kan man ha nytta av Mathematica för att kontrollera att man följer kokboken.
ode
y '' x
y x
2y x
y x
A
A
C x2
Bx
2 y' x
3yx
Bx
3y x
2x
2
2 x2
5
2 x2
5
5
C x2
ode
3 A
Bx
C x2 
2 B
2C x
2C
HH/ITE/BN
id
Ordinära differentialekvationer och Mathematica
CoefficientList
3A
2B
2 C, 3 B
,x &
17
ode
4 C, 3 C
5, 0, 2
Solve id
17
A
8
,B
2
,C
27

9
3
Exempel: Sök y p x då y ''
6 y'
5y
4x
9
är given.
Lösningsförslag: Först bestämmer vi den homogena lösningen med hjälp av karakteristiska ekvationen
r2
6r
Solver2
r
6r
1, r
5
5
0
6
2
r1,2

6 2

2
5
3
r1
r2
2
5
1
Fall 1
yh x
C1
5x
C2
x
0
5
Sedan en partikulärlösning.
Steg 1: Vi har f x
4x
9
Steg 2: Vi ser att y p x
så ansätt enligt tabell y p x
A
4x
.
yh x .
Steg 3: Identifiera koefficienter. Först bäddar vi upp med derivator
yp x
4x
A
yp ' x
4A
4x
y p '' x
16 A
4x
Insättning i differentialekvationen
16 A
4x
6 4A
4x

5 A
4x

9
4x
och identifiering av koefficienter framför varje funktionstyp (FT) i vänsterledet (VL) och högerledet (HL)
FT
4x
VL
16 A
HL
24 A
5A
9
Detta linjära ekvationssystem ska alltid ge en entydig lösning. Om inte så är ansatsen av y p x felaktig eller felräkning på vägen! Gå
i så fall tillbaka till Steg 1! Här får vi
Solve
16 A
24 A
5A
9
1
A

5
Så y p x
1
5
4x
varav allmänna lösningen som summan av den homogena lösningen och partikulärlösningen
yx
yh x
yp x
C1
5x
C2
x
1
5
4x
Med Mathematica har vi direkt
DSolvey '' x
6 y' x
5y x
9
4x
, y x , x
4x
y x
c1
x
c2
5x

5
Exempel: Sök y p x då y ''
4 y'
13 y
2 cos 4x
8 sin 4x är given.
Lösningsförslag: Först bestämmer vi den homogena lösningen med hjälp av karakteristiska ekvationen
r2
4r
13
0
r1,2
4
2
4 2
2
13
2
3
Fall 3
yh x
2x
C1 cos 3x
C2 sin 3x
18
Ordinära differentialekvationer och Mathematica
Solver2
r
2
4r
13
3 , r
2
HH/ITE/BN
0
3
Sedan en partikulärlösning.
Steg 1: Vi har f x
2 cos 4x
Steg 2: Vi ser att y p x
8 sin 4x så ansätt enligt tabell y p x
A cos 4x
B sin 4x .
yh x .
Steg 3: Identifiera koefficienter. Först bäddar vi upp med derivator
yp x
A cos 4x
B sin 4x
yp ' x
4 A sin 4x
y p '' x
16 A cos 4x
4 B cos 4x
4 A sin 4x
4 B cos 4x
16 B sin 4x
Insättning i differentialekvationen
16 A cos 4x
16 B sin 4x
4
13 A cos 4x
B sin 4x
2 cos 4x
8 sin 4x
och identifiering av koefficienter framför varje funktionstyp (FT) i vänsterledet (VL) och högerledet (HL)
FT
VL
cos 4x
sin 4x
16 A
16 B
HL
16 B
16 A
13 A
13 B
2
8
Detta linjära ekvationssystem ska alltid ge en entydig lösning. Om inte så är ansatsen av y p x felaktig eller felräkning på vägen! Gå
i så fall tillbaka till Steg 1! Här får vi
Solve
16 A
16 B
122
A
122
265
2,
16 B
16 A
13 B
8
56
,B

265
Så y p x
13 A
265
56
265
cos 4 x
sin 4 x varav allmänna lösningen som summan av den homogena lösningen och partikulärlösningen
yx
yh x
2x
yp x
C1 cos 3x
122
265
C2 sin 3x
cos 4x
56
265
sin 4x
Med Mathematica har vi direkt
DSolve y '' x
y x
c1
2x
4 y' x
sin 3 x
Exempel: Sök y p x då y ''
c2
6 y'
13 y x
2x
2 Cos 4 x
56
sin 4 x
cos 4 x 
265
5
3x
Simplify
122
cos 3 x
9y
8 Sin 4 x , y x , x
265
är given.
Lösningsförslag: Först bestämmer vi den homogena lösningen med hjälp av karakteristiska ekvationen
r2
Solver2
r
6r
3, r
6r
9
9
0
r1,2
6
2
6 2
2
9
2
0
Fall 2
yh x
C1 x
0
3
Sedan en partikulärlösning.
Steg 1: Vi har f x
5
Steg 2: Vi ser att y p x
3x
A
så ansätt enligt tabell y p x
3x
yh x
yp x
Ax
A
3x
3x
.
yh x
yp x
A x2
Steg 3: Identifiera koefficienter. Först bäddar vi upp med derivator av en produkt
3x
yh x .
C2
3x
HH/ITE/BN
Ordinära differentialekvationer och Mathematica
A x2
yp x
19
3x
3x
yp ' x
A 2 x
y p '' x
A 3 2x
x2
3x
3
A  3 x2
2
3x
 3 x2
3x
2 x
3x
2 x
A 9 x2
3
12 x
3x
2
Insättning i differentialekvationen
A 9 x2
12 x
3x
2
6 A  3 x2
2 x
3x
9 A x2
3x
5
3x
och identifiering av koefficienter framför varje funktionstyp (FT) i vänsterledet (VL) och högerledet (HL)
FT
VL
2
x
3x
x
3x
9A
18 A
12 A
3x
HL
9A
0
12 A
0
2A
5
Detta linjära ekvationssystem ska alltid ge en entydig lösning. Om inte så är ansatsen av y p x felaktig eller felräkning på vägen! Gå
i så fall tillbaka till Steg 1! När vi har korrigerat y p x genom multiplikation med x får vi lika många urartade ekvationer 0
0 som
antal korrigeringar. Bra kontroll! Här är det de två första.
Solve
9A
18 A
9A
0,
12 A
12 A
0, 2 A
5
5
A

2
5
2
Så y p x
x2
3x
varav allmänna lösningen som summan av den homogena lösningen och partikulärlösningen
yx
3x
Mathematica noterar naturligtvis att både
gen därefter.
DSolvey '' x
y x
c2
3x
6 y' x
x
3x
c1
5
3x
yp x
3x
och x
9y x
3x
5
C1 x
C2
3x
5
2
x2
3x
ingår i homogena lösningen och korrigerar ansatsen till partikulärlösnin-
, y x , x
x2 
2
y'' 2 y ' 3 y
y0 1
y' 0 0
Exempel: Lös BVP
yh x
x2
ODE
.
BV
Lösningsförslag: Först bestämmer vi den homogena lösningen med hjälp av karakteristiska ekvationen
r2
Solver2
r
3, r
2r
3
2r
0
3
r1,2
2 2
2
2
2
3
1
r1
r2
2
3
1
Fall 1
0
1
Sedan en partikulärlösning.
Steg 1: Vi har f x
x2 så ansätt enligt tabell y p x
Steg 2: Vi ser att y p x
A
Bx
C x2 .
yh x .
Steg 3: Identifiera koefficienter. Först bäddar vi upp med derivator
y p x A B x C x2
yp ' x B 2 C x
y p '' x
2C
Insättning i differentialekvationen
2C
2 B
2C x
3 A
Bx
C x2 
x2
yh x
C1
3x
C2
x
20
Ordinära differentialekvationer och Mathematica
HH/ITE/BN
och identifiering av koefficienter framför varje funktionstyp (FT) i vänsterledet (VL) och högerledet (HL)
FT
VL
2
x
HL
3C
1
4C
x
x0
2C
1
3B
2B
0
3A
0
Detta linjära ekvationssystem ska alltid ge en entydig lösning. Om inte så är ansatsen av y p x felaktig eller felräkning på vägen! Gå
i så fall tillbaka till Steg 1! Här får vi
Solve
3C
1, 4 C
3B
4
1
14
A
,B
,C
27
14
27
Så y p x
x
1
3
2B
3A
0

9
4
9
0, 2 C
3
x2 varav allmänna lösningen som summan av den homogena lösningen och partikulärlösningen
yx
yh x
yp x
3x
C1
14
27
x
C2
4
9
1
3
x
x2
y0 1
. Frestas inte att sätta in (BV) redan i yh x ! Det blir enklare
y' 0 0
Nu är det dax att bestämma C1 och C2 med hjälp av BV
räkningar men garanterat fel! (BV) gäller ju för lösningen till (ODE) inte dess homogena variant. Så först derivatan av y x
y' x
3x
3 C1
C2
4
9
x
2
3
x
sedan de efterlängtade
y0 1
:
y' 0 0
BV
30
C1
C2
30
3 C1
0
C2
14
27
0
4
9
4
9
1
3
0
2
3
0
02
0
1
C1
C2
3 C1
41
27
C2
C1
4
9
C2
29
108
5
4
Med Mathematica har vi direkt
DSolvey '' x
x2
4x
3
9
2 y' x
29
3x
x2 , y 0
3y x
x
5
1, y ' 0
0, y x , x
Expand
14
y x

108
4
27
Notera ännu en gång att begynnelsevärden (BV) och randvillkor (RV) skall appliceras på den allmänna lösningen y x inte på
den homogena lösningen yh x !
HH/ITE/BN
Ordinära differentialekvationer och Mathematica
21
ť Modellering av (ODE) för tillämpningar - några spridda tips och exempel
Att modellera verkliga problem med hjälp av matematik kräver erfarenhet och ett öppet sinnelag. Till sin natur är denna verksamhet
oftast helt skild från den matematikundervisning man möter i skolan och upplevs av naturliga skäl som svår. Inte minst på grund av
att man ska hämta kunskap och inspiration från flera grundläggande ämnen, inte bara matematik även om detta kommer att bli själva
språket. Det gäller att "se" på verkligheten genom glasögon som identifierar och tydliggör fenomen som kan kläs med matematiska
begrepp som derivata, integral och differentialekvation. Med detta betraktelsesätt tillägnar sig en modern ingenjör en klar konkurrensfördel eftersom modellering och simulering är ett absolut krav för att klara de allt kortare ledtider som krävs för att utveckla nya
produkter. Detta må vara allt från bilar, mobiltelefoner, medicinsk utrustning och läkemedel till schemaläggning av flygplan och
kommer att accentueras alltmer i framtiden. Se vidare i e-boken "Något om Matematisk modellering och Mathematica".
Vi börjar med några modelleringstips och avslutar med typiska fysikaliska fenomen som modelleras med (ODE). För fler exempel
och även mera avancerade sådana hänvisas till det senare avsnittet "Blandade exempel för resten". Den för ingenjörer så vanliga
Newtons accelerationslag behandlas utförligt med exempel i e-boken "Något om Mekanik-Dynamik och Mathematica".
Läsa spelet
Ofta gäller det att kunna dechiffrera en verbal beskrivning till en matematisk beskrivning. Vid differentialekvationer gäller det att
först bestämma vad som är den oberoende variabeln, i många praktiska tillämpningar är det tiden t, sedan den beroende variabeln,
det vill säga den sökta funktionen, säg y t .
Det viktiga är att formuleringen av en (ODE) ska göras och gälla för en godtycklig tidpunkt t i det intervall där man är intresserad
av att studera systemet, exempelvis t 0. Begynnelsevärden (BV) fixerar sedan alla Ci .
Som hjälp vid dechiffreringen är det bra att känna igen vanliga formuleringar som skall tolkas som derivata av den sökta funktionen
y
Ändring av y per tidsenhet
t
y
Ändringshastigheten av y
t
eller y' t .
eller y ' t .
När det sedan gäller att formulera själva (ODE):n är i en uppsjö av tillämpningar olika typer av proportionalitet inblandade. Vissa
ord, eller kombinationer av ord, kan också direkt associeras med matematiska krumelurer, exempelvis
y t är proportionell mot t
yt
kt.
k
y t är omvänt proportionell mot t
yt
k
.
t
k
Ändring av y t per tidsenhet är proportionell mot y t
y' t
y' t
ky t .
k
och t
Ändring av y t per tidsenhet är proportionell mot y t 
y' t
y' t
k y t t.
k
Ändring av y t per tidsenhet är proportionell mot y t 
och skillnaden mellan m och y t
y' t
y' t
ky t m
och omvänt proportionell mot skillnaden mellan m och y t
y t :s hastighet är proportionell mot y t 
y' t
k
yt .
m yt
k
k
y' t
ky t
m yt
.
m yt
Efter man löst (ODE) och bestämt alla Ci med hjälp av (BV) återstår avslutningsvis att fixera proportionalitetskonstanten k och
kanske andra konstanter med hjälp av några i problemtexten beskrivna ögonblicksbilder (RV) under resan.
Fysikaliska principer
Ofta kommer någon fysikalisk lag eller princip till användning, exempelvis
Newtons accelerationslag: m
2y
t2
F eller m y
F. Förmodligen den viktigaste ekvationen för en ingenjör. Se vidare
"Blandade exempel för resten" och speciellt i e-boken "Något om Mekanik-Dynamik med Mathematica" som
helt tillägnas denna ekvation. Vid handräkning brukar det ofta vara populärt att med kedjeregeln skriva om
andraderivatan y
def
y
y
y
t KR
t
y
y
y
y
för att få en första ordningens (ODE), inte sällan separabel.
22
Ordinära differentialekvationer och Mathematica
Energins bevarande
massans oförstörbarhet. Massa
HH/ITE/BN
volym densitet (eller koncentration).
Arkimedes princip: Då en kropp nedsänkes i en vätska påverkas denna av en lyftkraft som är lika stor som den
undanträngda vätskans tyngd.
SI-enheter
Håll koll på alla enheter. Räkna alltid i SI-enheter! Se vidare i e-boken "Något om SI-enheter och Mathematica".
Dimensionsanalys
Se till att det är lika många "äpplen" och "päron" på båda sidor om – tecknet i (ODE):n! Dimensionsanalys är ett utmärkt stöd vid
modellering och ett studium av de storheter som ingår i problemet ger nästan alltid direkta tips som leder till målet. Se vidare i eboken "Något om Dimensionsanalys och Mathematica".
Radioaktivt sönderfall, bakteriekultur, kapitalisering av pengar
Väldigt många problem beskrivs med den enkla proportionalitetsmodellen att
tillväxthastigheten vid varje tidpunkt är proportionell mot aktuell mängd. Detta
brukar ibland kallas Malthus lag.
Om vi som exempel tar radioaktivt sönderfall och låter aktuell mängd vid tiden t vara m t får vi direkt
m' t
m0
BVP
Detta är en separabel (ODE)
m
m
k t
1
m
m
k t
C1
Men den är också en linjär första ordningens (ODE) m ' t
t

kt
kt
mt
mAvt
m t
0
{Separabel}
DSolve m ' t
kt
c1
m t
DSolve
kt
M
mt
ln m
km t
0 t
km t , m t ,t
ODE
BV
kt
C1
k t C1
mt
mt
C1
0 med med integrerande faktorn IF
kt
mt
C1
mt
C1
kt
kt
mt
k t
k t,
C1
k t.
så
. Eller med Mathematica
First

Sedan bestäms C1 av (BV) m 0
mAvt

kt
km t
M
m' t
M: M
C1
k0
C1
km t , m 0
M . Eller hela resan med Mathematica
M ,m t ,t
First

Vid radioaktivt sönderfall brukar proportionalitetskonstanten k fixeras antingen av en ögonblicksbild, t.ex. m Τ
till halveringstiden H.
mΤ , eller kopplas
M
ekv
m t
. mAvt . t
H
2
M
M
Hk
2
kVärde
Solve ekv, k, Reals
First
log 2

k
H
Vi ser att k 0 som sig bör eftersom mängden minskar (sönderfaller). Man kunde naturligtvis lika gärna satt " k" i (ODE):n ovan. I
problemtexter brukar halveringstiden vara given och k bestäms som ovan, eller så är k given och halveringstiden söks. Proportionalitetskonstanten k avgör hur snabbt det radioaktiva ämnet klingar av. Ju rödare kurvan är i figuren nedan ju större är k och ju fortare
går det!
HH/ITE/BN
Ordinära differentialekvationer och Mathematica
kt
PlotEvaluateTable
PlotStyle
Hue
, k, 10 , t, 0, 5 , PlotRange
Range 0.1, 1, 0.1 , AxesLabel
23
All,
"t", "m t " 
mt
1.0
0.8
0.6
0.4
0.2
1
2
3
4
5
t
Fler exempel ges under det senare avsnittet "Blandade exempel för resten".
Newtons avsvalningslag
Kroppar, exempelvis kaffe, kaka från ugnen eller metallstycke, som ställs till avsvalning
eller för all del uppvärmning följer den enkla lagen att temperaturändringen per tidsenhet
är proportionell mot skillnaden mellan omgivande temperatur och aktuell temperatur
i kroppen. Omgivningen anses vara så stor att dess temperatur är konstant under
resans gång. Detta samband kallas Newtons avsvalningslag.
Lösningsförslag: Om vi antar att temperaturen i kroppen är T t vid tiden t får vi direkt
BVP
T
Tomg T
Detta är en separabel (ODE);
Tomg
T
k t C1
T 0
T0
T t
T
k t
C1
Men också en linjär första ordningens (ODE); T ' t
kT t
kTomg med IF
kt
T t
kt
kTomg t
kt
T t
1
kTomg  k
Slutligen återstår bara att fixera C1 med hjälp av (BV) T 0
C1
ln Tomg
t
kt
C1
k t.
Tomg

ODE
BV
1
Tomg T
k t, som har lösningen
k t C1
Tomg
k Tomg
T t
Separabel
T t
T' t
kt
T0 : T0
allmänna lösningen T t ovan så har vi lösningen till (BVP) T t
C1 
k t
T t
Tomg
Tomg
Tomg
Tomg 
så
C1
kt
t
C1
k0
C1
T0
k t,

kt
T t
kt
.
T0
kt
kTomg
Tomg , och att meka in den i
.
Typiskt för alla lösningskurvor är att de alla startar i T0 och närmar sig asymptotiskt Tomg ju längre tiden lider. Gäller både avsvalning och uppvärmning! Detta inses av (ODE) med k
0, ty Tomg
så Tomg är en asymptot. På analogt sätt Tomg
T' t
T t
T t
T' t
0 ger avsvalning men då T t
Avsvalning av kaffe i rumstemeratur eller uppvärmning av drink på stranden.
Tk
T t
Td
T t
DSolve
60
kt
DSolve
20
13
T' t
T' t
T 0
k 25
85
T t
ODE
BV
T' t
T 0
respektive
k 25
T t
,T 0
85 , T t , t
k 20
T t
,T 0
7 ,T t ,t
k 20
7
T t
ODE
BV
Simplify
25
T' t
kt
T' t
0 har vi uppvärmning mot Tomg . Som vanligt avgör k hur snabbt
processen förlöper.
BVP
Tomg
Simplify

Ju rödare kurvan är i figuren nedan ju större är k och ju fortare går processerna!
Plot Evaluate Table T t
. Tk, Td , k, 10
, t, 0, 1 , PlotRange
0, 100 ,
PlotStyle Hue
Range 0.1, 1, 0.1 , AxesLabel
"t", "T t
C "
0
24
Ordinära differentialekvationer och Mathematica
T t
100
HH/ITE/BN
C
80
60
40
20
0.0
0.2
0.4
0.6
0.8
1.0
t
Fler exempel ges under det senare avsnittet "Blandade exempel för resten".
Torricellis lag
Hastigheten som en vätska strömmar ut genom ett hål i en tank är v
2gh ,
där g är tyngdaccelerationen och h vattendjupet vid hålet. Detta sammband
kallas Torricellis lag. Så om hålets area är A har vi flödet q Av. I praktiken
har vi energiförluster vid hålet så det teoretiska uttrycket för v måste korrigeras
med en faktor Η 1 som beror på hålets utseende. För ett cirkulärt hål är
Η 0.6 om det strömmar rakt ut i luften.
Så med ord har vi Torricellis lag: Då en vätska strömmar med flödet q ur en tank gäller att detta i varje ögonblick är proportionellt
mot kvadratroten ur vattendjupet h. Eftersom flöde är volymändring per tidsenhet kan vi sammanfatta i en (ODE)
V
t
q
k h
Där V V t är aktuell vätskevolym. Vattendjupet h t ska betraktas som en koordinat för vattenytans läge, enligt figur. Lagen
brukar ibland formuleras med minustecken framför proportionalitetskonstanten k, för att redan vid modelleringen återspegla det
faktum att volymen minskar i tanken. Speciellt viktigt är det naturligtvis att det blir rätt tecken om k är given numeriskt. Problemet är
i allmänhet att för givet utseende på tanken beskriva V
V
t
V h varav sedan
V
h
h
t
med kedjeregeln eller implicit derivering ger
en (ODE) som bestämmer h t . Proportionalitetskonstanten k fixeras ofta med en ögonblicksbild någon gång under tömningen.
Exempel: En tank i form av en stående cylinder är helt fylld med
vatten. Så öppnas en kran i botten så att vattnet strömmar ut med
en hastighet som i varje ögonblick är proportionell mot kvadratroten
ur vattendjupet, Torricellis lag. Hur lång tid tar det att tömma
tanken om den efter T s är tömd till hälften?
Lösningsförslag: Om R är tankens radie och h en koordinat för vattendjupet räknat från botten har vi
säga ΠR2 h ' t
k ht
h' t
k
ΠR2
ht
h' t
t
ΠR2 h
k h . Det vill
k h t , med sedvanlig omdöpning av konstanter. Antag att tanken har
höjden H så kan vi sammanfatta
BVP
h' t
h0
hT
k ht
H
ODE
BV
H
2
RV
Vi börjar med handräkningsversionen och obestämda intergraler för en separabel (ODE). Den är inte linjär och vi börjar som vanligt
med allmänna lösningen
h' t
k ht
h
12
h
k t
C1
Nu återstår att bestämma konstanten C1 med hjälp av (BV) h 0
1
12 1
H: H
h
12 1
1
4
k 0
kt
C1
C1
2
1
kt
4
ht
C1
2
2
C1 
H . Sedan fixerar vi k med hjälp
HH/ITE/BN
Ordinära differentialekvationer och Mathematica
H H
:
2 2
av (RV) h T
1
4
k T
2
2
H 
H
T
k
 2
H
2 . Så till slut h t
H
Tömningstiden Τ får vi sedan genom att lösa ekvationen 0
4 T2
 2
4 T2
 2
2
2 Τ
25
2 T
Τ
2 t
2
2
2 T .
2  T. Avslutningsvis låter vi
Mathematica lösa det ekvationssystem som formas av separation och bestämda integraler med (BV) och (RV) som gränser
1
H 2
ekv
H

T
2
0
h
2
H
k T, 2
2
Τ
h
H
H
Solve ekv, k, Τ
 2
1
0
 k t, 
h

 k t
0
h
k Τ
Simplify
H
k
,Τ
2
2  T
T
Till slut den moderna vägen med Mathematica. Som vanligt är den oberoende variabeln, här tiden t, den sammanhållande länken i
en (ODE). Vi har en första ordningens (ODE), så vi kan bara ta med ett av (BV) eller (RV), man brukar välja (BV).
hAvt
DSolveh ' t
k
1
h t
H kt
4
4H
h t
,h 0
H, h t , t
Last
k 2 t2 
4
Sedan k med (RV) och tömningstiden Τ
H
 , 0 . hAvt . t
2
Solveh t
 2
2
T, Τ , k, Τ 
FullSimplify
First
H
,Τ
k
2
2  T
T
Typiska arbetsgången med DSolve följt av Solve. Tillämpad matemaik är typiskt att formulera och lösa (system av) ekvationer.
Ytterligare exempel ges under det senare avsnittet "Blandade exempel för resten".
Blandningsproblem
I figuren åskådliggörs ett kontinuerligt blandningsförlopp. Det finns gott
om tillämpningar som kan modelleras på detta sätt, exempelvis förorening
i en tank, ett organ som genomströmmas av blod, en sjö eller våtmark med
anslutande åar. Uppgiften är att studera hur koncentrationen c t av något
ämne i tanken utvecklas över tiden. I figuren är följande definierat. Vi
räknar här i SI–enheter. Bra vana, men det viktiga är att enheterna som
används nedan är konsistenta
V t
Volym i tanken, m3 .
cin t
Koncentration av ämnet i qin t , kg m3 .
ct
Koncentration av ämnet i tanken, kg m3 .
qut t
Volymflödet ut ur tanken, m3 s .
qin t
Volymflödet in i tanken, m3 s .
cut t
Koncentration av ämnet i qut t , kg m3 .
Vi tar hjälp av den fysikaliska principen om massans oförstörbarhet och studerar tanken under en liten tid t vid godtycklig tidpunkt
t. Då måste det gälla
ökning av massa i tanken
massa in
massa ut
Nu gäller det att översätta detta till matematik innehållande storheterna i figuren. Här kommer dimensionsanalys väl till pass.
ökning av massa i tanken
V t c t
massa in
qin t cin t t
massa ut
qut t cut t t
m3  kg m3 
m3 s kg m3  s
m3 s kg m3  s
kg
kg
kg
26
Ordinära differentialekvationer och Mathematica
HH/ITE/BN
Ok, vi har [kg] på båda sidor! Dividera båda sidor med t, gå i gräns och känn igen derivatans definition.
V t ct
t
qin t cin t
qut t cut t
V t ct
t
t 0
qin t cin t
qut t cut t
Nu brukar man förutsätta perfekt omrörning. Detta innebär att cut t c t . Man kan se det som att man upprepat "släpper i en liten
droppe, rör om ordentligt och släpper ut en liten droppe". Vi sammanfattar den färdiga formuleringen av blandningsproblemet
V t ct
t
BVP
c0
där V t
qin t cin t
qut t c t
ODE
c0
BV
t
q
0 in
V 0
t
qut t
t
Hur svår denna (ODE) blir att lösa beror på de funktioner som ingår i högerledet. Inte sällan är cin t cin konstant under den tid vi
ämnar studera systemet. Detsamma brukar gälla för in- och utflödena, annars skulle volymen i tanken öka eller minska för att så
småningom bli helt tom. Jämför en sjö eller ett organ i kroppen.
Exempel: En sjö har volymen 105 m3 . Vid en tidpunkt uppmättes
koncentrationen kvicksilver i sjön till 10 mg m3 . Från en å rinner det
in vatten med flödet 2 m3 h innehållande kvicksilver med 3 mg m3 .
Sök koncentrationen kvicksilver i sjön som funktion av tiden efter
uppmätningen. Anta perfekt omrörning samt att det finns ett utlopp
från sjön så att dess volym är konstant över tiden.
Lösningsförslag: Typiskt blandningsproblem. Med kända beteckningar har vi
t
BVP
V t c t 
c 0
qin t cin t
qut t c t
BV
c0
105
V t konstant 105 , qin qut 2
ODE
cin 3, perfekt omrörning
c
t
2 3
c 0
cut c
2 c
10
ODE
BV
Denna kan lösas som linjär eller separabel. Vi börjar känna igen att det är "samma" (ODE) antingen vi har blandningsproblem,
Newtons avsvalningslag, radioaktivt sönderfall Matematik är användbart!
Vi provar med linjär 105
c
t
6
t

2c
5
2 10
t
c' t
2 10
2 10
ct
2 10
5
t
5
t
5
6 10
ct
6 10 5 . Integrerande faktorn IF
ct
5
6 10
5
2 10
5
2 10
Förväxla inte c t med C1 Den senare fixeras av (BV) c 0
ct
3
2 10
7
cAvt
c t
5
t.
5
t
10: 10
5
t
C1
3
2 10
ct
ct
C1
5
3
C1
2 10
5
0
t
6 10
2 10
C1
5
5
5
t
C1
t
2 10
5
t
. Så
t
7. Så slutligen lösningen till (BVP)
Med Mathematica har vi direkt
DSolve105 c ' t
7
2 10

2 10
t 50 000
2 3
c t
,c 0
10, c t , t
Simplify
3
Å så här ser skådespelet ut under de trettio första åren
PlotEvaluate
3, c t
PlotStyle
. cAvt . t
Red, Dashing
0.02
c t mg m3
10
8
6
4
2
0
5
10
15
20
25
30
t år
24
365 t , t, 0, 30 , PlotRange
, Orange , AxesLabel
"t
0, 10 ,
år ", "c t
mg m3 "
HH/ITE/BN
Ordinära differentialekvationer och Mathematica
27
Vi ser både i den analytiska lösningen och i figuren ovan att föroreningen i sjön efter lång tid närmar sig den i inloppet. Verkar ju
rimligt. Så är det alltid !
Exempel: Samma förutsättningar som föregående exempel med enda skillnaden att det rinner in rent vatten så att sjön tvättas ren.
Hur länge dröjer det innan koncentrationen har sjunkit till hälften?
Lösningsförslag: Ny cin t
0, så vi har att ta ställning till
105
BVP
c
t
c0
2 0
2 c
ODE
10
BV
c
t
Denna kan lösas som linjär eller separabel. Vi provar med separation, 105
ln c
2 10
5
t
C1
lösningen till (BVP) c t
cAvt
c t
ct
2 10
10
2 10
C1
5
t
t 50 000
t
. Sedan fixeras C1 av (BV) c 0
10: 10
C1
c
2 10
5
5
2 10
0
C1
t
C1
10, varav slutligen
. Med Mathematica har vi direkt
DSolve105 c ' t
10
5
1
c
2c
2
0
2c t , c 0
10, c t , t

Halveringstiden i timmar.
10
Solvec t
. cAvt, t, Reals
2
t
50 000 log 2
Å så här ser skådespelet ut under de trettio första åren
PlotEvaluate c t
PlotRange
. cAvt . t
All, AxesLabel
24
365 t , t, 0, 30 , PlotStyle
"t
år ", "c t
Orange,
mg m3 "
c t mg m3
10
8
6
4
2
5
10
15
20
25
30
t år
Exempel: Avslutningsvis ursprungsexemplet med enda skillnaden att
en bäver har mekat om sitt bo i ingående å så att inflödet ökat till 6 m3 h
Lösningsförslag: Volymen i sjön blir nu inte konstant, utan V t
BVP
t
105
c0
t
q
0 in
V 0
4 t c t 
6 3
t
qut t
2 ct
10
t
105
6
2 t. Så
ODE
BV
Denna kan lösas antingen som linjär eller separabel. Vi provar med linjär. Först utvecklar vi vänsterledet som är derivatan av en
produkt, sedan stuvar vi om till standardform
t
105
4 t c t 
6
Integrerande faktorn IF

105 4 t
t
18
6
4
2c t
4c t
105
ln105 4 t
105
4 t . Så
32
4 t c' t
18
2c t
c' t
6
105 4 t
ct
18
105 4 t
28
Ordinära differentialekvationer och Mathematica
t
32
105
4 t
32
105
Fixera sedan C1 med (BV) c 0
32
105
ct
4 t
10: 10
4 t
18
12 1
ct
3
32
4 t
12 1 1
4 t
4
4 0
7105 
32
105
105
C1 105
3
ct
18
105 4 t
32
105
18 105
ct
ct
C1 105
3
12
4 t
4 t
t
C1
32
32
7105  . Så slutligen
C1
4 t
C1
HH/ITE/BN
32
ct
7
3
1 4 10
5
t
32
Vi låter Mathematica visa musklerna genom att direkt bestämma både c t och V t ur ett system av differentialekvationer.
cÅV
DSolveD V t c t , t
V' t
6
c 0
10, V 0
c t ,V t
4 t
3
2c t ,
105 ,
, t
8 750 000
V t
6
2,
25 000 , c t
t
25 000
ODE
blandning
ODE
volym
BV
Simplify
10
3
32
Å så här ser skådespelet ut under de trettio första åren. Vi ser i lösningen ovan och figuren nedan att c t
cin
3, som sig bör.
V t
PlotEvaluatec t ,
 . cÅV . t
24
365 t, t, 0, 30 , PlotRange
0, 12 ,
105
PlotStyle
Orange, Blue , AxesLabel
"t
mg m3 ,V t
år ", "c t
105 m3 "
c t mg m3 ,V t 10 5 m3
12
10
8
6
4
2
0
5
10
15
20
25
30
t år
Ytterligare exempel ges under det senare avsnittet "Blandade exempel för resten".
ť Numerisk lösning av (ODE)
Vi ska använda en så kallad Finita differensmetod (FDM) för att lösa en differentialekvation. Detta är en metod för att söka en
approximativ numerisk lösning till ett begynnelse- eller randvärdesproblem som av olika skäl är så besvärligt att vi inte klarar av att
lösa det analytiskt med metoder vi känner från tidigare moment i kursen. Detta är den "vanligaste" situationen för en ingenjör i
y ' 3 y 1 ODE
verkliga livet! Vi testar den på ett enkelt (BVP)
, så har vi något att jämföra med.
y0 0
BV
yAvx
DSolve
y' x
3y x
1, y 0
0 ,y x ,x
First
Simplify
3x
1
y x

3
3
Plot y x
. yAvx, x, 0, 1 , PlotStyle
y
0.30
0.25
0.20
0.15
0.10
0.05
0.2
0.4
0.6
0.8
1.0
x
Blue, AxesLabel
"x", "y"
HH/ITE/BN
Ordinära differentialekvationer och Mathematica
29
Vi söker nu en approximativ numerisk lösning i ett område, säg x 0, 1 . Idén är att dela in området i ett valfritt antal mindre
intervall, så kallade element. Gränspunkterna mellan elementen kallas noder. I de inre noderna ersätts sedan y ' med sin differenskvot. Vi får då ett ekvationssystem ur vilket y kan bestämmas i varje nod. Vi har därmed fått en approximativ numerisk lösning
som är styckvis linjär över elementen. Det typiska är nu att ju finare indelning vi gör ju bättre blir den approximativa numeriska
lösningen. Metoden kan förfinas till den så kallade Finita elementmetoden (FEM) vilken är det helt dominerande redskapet för
datorbaserad beräkning av hållfasthet för en maskindetalj, lyftkraften hos en flygplansvinge eller väderprognoser.
Antag att vi delar in vårt område i n st lika långa element. Elementen behöver inte vara lika långa, snarare tvärtom, där det "händer"
mycket måste man använda mindre element för att kunna följa lösningen. Omvänt gäller då större element i "lugnare" områden. Vi
börjar med en grov indelning som har stora element. Dessa blir då h långa.
1.0
n
3; h
n
0.333333
Vi har alltså indelningen med nodernas lägen xi och de eftersökta nodvärdena yi . Här vet vi på grund av (BV) att y 0
x0
x0
0
x1
xi
xn
x
1
x1
x2
0
x0
x1
xn
y0
0.
y0 0
h
y1
h
y2
1
yn
Approximera nu derivatan i varje inre nod med den så kallade bakåtdifferenskvoten, det vill säga yi '
yi yi
1
h
. Detta leder till ett
ekvationssystem där varje rad motsvarar differentialekvationen applicerad i tur och ordning på de inre noderna.
yi
ekv
yi
1
Table
3 yi
1, i, 1, n 
h
3 y1
3. y1
y0
1, 3 y2
3. y2
y1
1, 3 y3
3. y3
y2
1
Nu är detta ekvationssystem singulärt och kräver liksom tidigare ett begynnelsevärde. Jämför bestämning av C1 för att fixera en
fysikalisk lösning till differentialekvationen! I vårt fall y 0 y0 0. Sedan är det bara att lösa det.
lösning
y1
Solve ekv . y0
0.166667, y2
0.25, y3
0
First
0.291667
eller i tabellform xi , yi .
i
Table , yi , i, 0, n  . lösning . y0
n
xyFDM
0
1
3
2
3
0
0
0.166667
0.25
1 0.291667
Slutligen en jämförelse med den analytiska lösningen.
Plot y x
. yAvx, x, 0, 1 , PlotStyle Blue, PlotRange All, AxesLabel
Epilog
Red, PointSize 0.025 , Line xyFDM , Point xyFDM
y, yFDM
0.30
0.25
0.20
0.15
0.10
0.05
0.2
0.4
0.6
0.8
1.0
x
"x", "y, yFDM" ,
30
Ordinära differentialekvationer och Mathematica
HH/ITE/BN
Som tidigare nämnts kan man förvänta sig bättre noggrannhet om en finare indelning väljs. I figuren nedan ses en jämförelse mellan
analytisk lösning och en indelning med 25 element. Man kan visa att felet i varje nod är proportionellt mot h.
y, yFDM
0.30
0.25
0.20
0.15
0.10
0.05
0.2
0.4
0.6
0.8
1.0
x
Metoden ovan är en variant av som brukar kallas Eulers stegmetod eftersom man egentligen inte behöver lösa ett ekvationssystem.
Metoden kan användas då derivatan kan räknas ut explicit y ' f x, y . I vårt modellexempel ovan är tydligen detta fallet, eftersom
y'
1
3 y. Med framåtdifferens får vi
yi
1
yi
1
h
3 yi
yi
1
yi
h1
3 yi , vilket visar att varje ny punkt på kurvan kan
xi , yi
h 1, 1 3 yi . Eller mer allmänt
beräknas genom att stega fram i tangentens riktning med start i (BV); xi 1 , yi 1
xi , yi
h 1, f xi , yi . Jämfört med FDM hamnar vi på andra sidan om den analytiska lösningen
xi 1 , yi 1
1.0
n
25; h
i
; xyEuler
1
Table
n
, If i
1, y0
0, yi
1
yi
h 1
3 yi
, i,
1, n
1 ;
n
Plot y x
. yAvx, x, 0, 1 , PlotStyle Blue, PlotRange All, AxesLabel
Epilog
Red, PointSize 0.025 , Line xyEuler , Point xyEuler
"x", "y, yEuler" ,
y, yEuler
0.30
0.25
0.20
0.15
0.10
0.05
0.2
0.4
0.6
0.8
1.0
x
Eulers metod kan trimmas till bättre noggrannhet genom att använda fler punkter när derivatan ska uppskattas. En sådan populär är
Runge-Kuttas metod. Vid högre ordningens differentialekvationer eller mer förfinad approximation av derivatan får vi emellertid
en starkare koppling mellan ekvationerna som då kräver ekvationslösning, t.ex. Gauss eliminationsmetod. Vi tar som exempel ett
andra ordningens randvärdesproblem
y '' 2 y '
y0 1
y1 2
RVP
y
1
x2
ODE
RV
med analytisk lösning
yAvx
DSolvey '' x
y 0
y x
x2
Plot y x
12
x 1
x
2 y' x
1, y 1
4x
6
x
x
y x
1
2, y x , x
1
. yAvx, x, 0, 1 , PlotStyle
2.0
1.8
1.6
1.4
1.2
0.4
0.6
0.8
First
Simplify
5
y
0.2
x2 ,
1.0
x
Blue, AxesLabel
"x", "y"
HH/ITE/BN
Ordinära differentialekvationer och Mathematica
31
Antag att vi delar in vårt område i n st lika långa element.
n
5;
Dessa blir då h långa.
1.0
h
n
0.2
Andraderivatan approximeras med framåtdifferenskvot och bakåtdifferenskvot yi ''
yi 1 yi
h
yi yi 1
h
yi
1
2 yi yi
h2
h
1
. Vi ser att denna blir
symmetrisk. Här har vi även möjlighet att göra approximationen av derivatan symmetrisk genom att använda centraldifferenskvot,
det vill säga yi '
yi
yi
1
1
2h
. Detta leder till ett ekvationssystem där varje rad motsvarar differentialekvationen applicerad i tur och
ordning på de inre noderna.
yi
ekv
1
2 yi
yi
yi
1
Table
1
yi
1
2
yi
h2
1
ih
2
, i, 1, n
y1 5. y2 y0 25. y0 2 y1 y2
y3 5. y4 y2 25. y2 2 y3 y4
0.96, y2 5. y3 y1 25. y1 2 y2 y3
0.64, y4 5. y5 y3 25. y3 2 y4 y5
Nu är detta ekvationssystem singulärt och kräver liksom tidigare randvärde (RV), här
lösning
y1
Solve ekv . y0
1.1, y2
1.24299, y3
1 
2h
1, yn
1.43733, y4
2
0.84,
0.36
y0
y1
y0
yn
1
. Nu är det bara att lösa det.
2
First
1.68898
eller i tabellform xi , yi .
i
Table , yi , i, 0, n  . lösning . y0
n
xyFDM
0
1
5
2
5
3
5
4
5
1
1, yn
2
1
1.1
1.24299
1.43733
1.68898
2
Slutligen en jämförelse med den analytiska lösningen.
Plot y x
. yAvx, x, 0, 1 , PlotStyle Blue, PlotRange All,
PlotRange
Automatic, 0, 0.4 , AxesLabel
"x", "y, yFDM" ,
Epilog
Red, PointSize 0.025 , Line xyFDM , Point xyFDM
y, yFDM
2.0
1.8
1.6
1.4
1.2
0.2
0.4
0.6
0.8
1.0
x
Som tidigare nämnts kan man förvänta sig bättre noggrannhet om en finare indelning väljs. Eftersom vi har symmetriska approximationer kan man visa att felet i varje nod är proportionellt mot h2 .
32
Ordinära differentialekvationer och Mathematica
HH/ITE/BN
ť Blandade exempel för resten
Exempel: En iskub som glömts på stranden smälter
så att volymändringen per tidsenhet är proportionell
mot dess area. Antag att sidan var 3 cm från början
och att den smält till 2 cm på 5 min. Hur länge dröjer
det innan den har smält bort ?
V
t
Lösningsförslag: Enligt uppgift har vi att volymändringen per tidsenhet är proportionell mot arean,
s3 och arean A
volymen V
s
t
3 s2
s
t
k 6 s2
6 s2 så med den populära kedjeregeln får vi
k. Direkt integration ger s
3
2
k 0
k 5
m
m
1
5
k
och m
kt
V
t
kA
KR
V
s
s
t
kA. Med sidan s blir
k 6 s2
s
s3 
s
t
k 6 s2
m och sedan med (BV) och (RV)
1
5
3, varav svaret på den brännande frågan: 0
t
3
t
15 min.
Avslutningsvis ett litet trick Om man vill få med randvillkoret s 5 2 kan man lösa ett system av (ODE) med vetskapen att k är
en konstant, det vill säga k ' t 0. Nu har vi totalt två ' så vi har möjlighet att ta med två BV
RV . Detta ger direkt s t och k.
sÅk
DSolveDs t
1
15
k t
3
, t
k t 6s t
2
, k' t
3, s 5
2, s t , k t
, t
t
,s t

10
5
med samma svar s t som ovan. Dock har vi "ursprungliga" k i
Solves t
t
0, s 0
V
t
kA. Smidigt! Slutligen smälttiden som ovan
0 . sÅk
15
Exempel: År 1988 kontaktade Vatikanen British Museum för att få hjälp med
att datera den så kallade Svepningen från Turin som påstods vara den svepning
som användes vid begravningen av Jesus från Nasaret. Den upptäcktes år 1356
och ansågs av många allt sedan dess vara äkta. British Museum bestämde sig för
att använde kol 14 metoden för att göra denna datering. Efter en undersökning
av fibrerna i svepningen visade det sig att mellan 92 och 93 av ursprungliga
mängden kol 14 fanns kvar. Kol 14 har halveringstiden 5730 år. Gör om den
analys som British Museum gjorde
Lösningsförslag: Om mängden kol 14 vid tiden t är m t så gäller för radioaktivt sönderfall att sönderfallshastigheten
m' t
k m t . Denna är separabel
mt
kt
C1
mAvt
m t
m
m
k t, varav
1
m
m
k t
.
DSolve m ' t
kt
c1
km t ,m t ,t
First

Randvillkoret med halveringstiden fixerar proportionalitetskonstanten k.
1
ekv
m t
c1
. mAvt . t
5730
2
c1
5730 k
c1
2
kVärde
Solve ekv, k, Reals
log 2
k

5730
First
C1
ln m
kt
C1
m
k t C1
m
C1
kt
HH/ITE/BN
Ordinära differentialekvationer och Mathematica
33
Lös ut tiden motsvarande en godtycklig andel p.
Simplify Solve m t
Τ
p c1
. mAvt . kVärde, t, Reals , p
0
First
5730 log p
t

log 2
Slutligen den ålder som motsvarar 92% och 93% av kvarvarande mängd jämfört med från början.
Τ .p
t
0.92, 0.93
689.286, 599.916
Svepningen kan därmed dateras till att vara från tidsepoken
1988
t .Τ .p
0.92, 0.93
Round
1299, 1388
Om man tror på kol 14 metoden kan alltså svepningen inte ha tillhört Jesus från Nasaret.
Exempel: En bakteriekultur dubbleras på 10 min. Hur lång
tid tar det innan kulturen har fem gånger så många bakterier
som från början ? Antag att tillväxthastigheten vid varje
tidpunkt är proportionell mot antalet bakterier. Detta
brukar ibland kallas Malthus lag.
Lösningsförslag: Efter översättning av sista meningen har vi modellen b ' t
1
b
b
k t
C1
bAvt
ln b
kt
DSolve b ' t
b t
kt
c1
C1
b
k t C1
kb t , b t ,t
b
C1 k t
bt
k b t . Denna är separabel,
C1
First

Kravet på dubblering efter 10 s fixerar proportionalitetskonstanten k.
ekv
c1
b t
10 k
2 c1
. bAvt . t
10
2 c1
kVärde
Solve ekv, k, Reals
First
log 2
k

10
Slutligen tiden till fem gånger så många bakterier som från början.
ekv
c1 2
b t
t 10
5 c1
. bAvt . kVärde
5 c1
Solve ekv, t, Reals
10 log 5
t

log 2
Exempel: En nackdel med tillväxten av bakteriekulturen i föregående exempel
är att antalet bakterier växer över alla gränser då t
. Detta kan korrigeras
med att bestämma en maximalt tillåten population som motsvarar mängden
tillgänglig mat och andra livsbegränsande betingelser. För att uppnå detta låter
vi högerledet också vara proportionellt mot skillnaden mellan maximala antalet
individer och aktuell mängd. Denna formulering brukar kallas Logistiska lagen
och vi exemplifierar med några andra kompisar.
kt
.
b
b
k t, med lösningen
34
Ordinära differentialekvationer och Mathematica
HH/ITE/BN
Lösningsförslag: Efter översättning av näst sista meningen har vi kaninmodellen k ' t Αk t K k t . Som tidigare indikerar Α
hur snabbt förökningen sker och K det maximala antalet kaniner som kan finnas i skogen, den så kallade kapaciteten. Denna (ODE)
är separabel,
kAvt
k
kK k
Α t, men kräver partialbråksuppdelning för att kunna integreras. Så vi låter Mathematica göra jobbet.
DSolve
k t
2
k' t
2K
ΑKt
ΑKt
K
Αk t
K
k t
,k 0
2 ,k t ,t
First

2
Ju rödare kurvan är i figuren nedan ju större är Α och ju fortare ökar populationen upp till fullbelagt, exempelvis 1000 st!
Plot Evaluate Table k t
. kAvt . K
PlotRange All, PlotStyle Hue
1000, Α, 0.01, 0.1, 0.01
, t, 0, 1 ,
Range 0.1, 1, 0.1 , AxesLabel
"t", "k t "
kt
1000
800
600
400
200
0.2
0.4
0.6
0.8
1.0
t
Exempel: En vanlig utveckling av logistiska lagen är med så kallad Allee effekt (Warder Clyde Allee 1885-1955, amerikansk
teoretisk ekolog) som tar hänsyn till att vissa arter har stora problem med att återhämta sig om populationen sjunker under en viss
nivå. På grund av genetisk degenerering kommer den då tyvärr att gå ett tragiskt öde till mötes och faktiskt dö ut. En art som tydligt
visade på denna effekt var vandringsduvan (Ectopistes migratorius), som efter passing by fick namnet passenger pigeon på
engelska. Under tidiga 1800-talet var den förmodligen jordens vanligaste fågel, och speciellt i nordamerika har det uppskattats att
det fanns enorma mängder, kanske 5-10 miljarder! När de europeiska emigranterna anlände gick det ryckte om flockar som var så
stora som 500 1.5 km, vilka det tog en hel dag för att passera ranchen! Med för kontinenten nya jaktvapen började så en förödande
jakt på fåglarna. Nyttojakten urartade till rena slaktorgier, se figur nedan till höger. Då förstår man att det knappt finns några
fotografier på den vackra fågeln, utan endast teckningar, som den till vänster nedan.
Effekten blev katastrofal för vandringsduvan som minskade kraftigt i antal mot slutet av 1800-talet. Massiva räddningsinsatser
visade sig verkningslösa, så den sista vandringsduvan, kallad Martha, dog den första september 1914 i Cincinatti Zoo.
Allee effekt monteras in i logistiska lagen som en extra faktor i högerledet som tar hand om den undre gränsen för överlevnad L. Om
vi startar ovanför denna gräns ökar populationen upp mot kapaciteten K och dör ut om vi startar under den. Om v t är antalet
vandringsduvor vid tiden t får vi alltså modellen v' t Αv t K v t v t L . Vi exemplifierar med ett normerat fall där K 1
och L 0.2. Vi ser tydligt Allee effekten när vi startar från lite olika begynnelsevärden på populationen.
pAvt
NDSolve
p' t
3p t 1 p t
p t
0.2 , p 0
,
p t , t, 0, 10
&
0.03 Range 15, 1, 1 ;
Plot Evaluate p t
. pAvt , t, 0, 10 , PlotRange All, AxesLabel
vt
1.0
0.8
0.6
0.4
0.2
2
4
6
8
10
t
"t", "v t "
HH/ITE/BN
Ordinära differentialekvationer och Mathematica
35
Exempel: Tomtens gröt har temperaturen 90 C då den lyfts från
spisen. Temperaturen sjunker sedan till 75 C på 5 min då temperaturen
i stugan är 20 C. Tomtens känsliga mage tål inget varmare än 50 C. Hur
länge måste han vänta ? Vad är temperaturen efter 15 min? Antag att
T
t
Newtons avsvalningslag
k Tomg
T gäller.
Lösningsförslag: Detta är en linjär eller separabel (ODE).
T t
Tomg
kt
C1
. Vidare har vi att Tomg
TAvt
DSolve
T t
kt
70
ODE
T' t
T t
Solve T t
1
k
k t
C1
ln Tomg
t
kt
C1
,T 0
20
20
k0
C1
70
kt
C1
70 C
. Samma svar med Mathematica
90 , T t , t
Simplify
First
20
Slutligen fixerar konstanten k av RV T 5
kVärde
90 : 90
BV blir T t
k 20
T
20 C. Först bestäms C1 av
(BV) T 0
Detta ger att lösningen till BVP
1
Tomg T
75 : 75
75 . TAvt . t
20
70
k5
k5
70
5, k, Reals
55
First
1
5
k
55
ln 70 
1
5
k
14
ln 11 .
Simplify
14

log
5
11
t
Slutligen svaret på Tomtens brännande fråga T t
Solve T t
log
50 . TAvt . kVärde
16 807
t
243
14
50 : 50
20
70
5
ln
14
11
7

30
ln 70 
t
5
14
ln 11 
t
5
ln 
ln
3
14
11
.

Simplify


log 
11
Direkt insättning i modellen ger svaret på sista frågan
T t
. TAvt . kVärde . t
15.
53.9541
Avslutningsvis ett litet trick Om man vill få med randvillkoret T 5 75 C kan man lösa ett system av (ODE) med vetskapen att k
är en konstant, det vill säga k ' t 0. Nu har vi totalt två ' så vi har möjlighet att ta med två BV
RV . Detta ger direkt T t och k,
det vill säga lösningen till ODE
BV
RV .
DSolve
T' t
T 0
1
k t
k t 20
90, T 5
14
,T t
log
5
T t , k' t
0,
75 , T t , k t , t
11t 5 14
5
1
Simplify
First
t
5
20
11
med samma svar som ovan. Newtons avsvalningslag kan lika gärna kallas för Newtons uppvärmningslag, ty T t
Tomg oberoende
av (BV). Hur snabbt detta går avgörs av k.
Plot Evaluate T t
. TAvt . k Range 0.1, 1, 0.1 , t, 0, 30 , PlotRange
0, 100 ,
PlotStyle Hue
Range 0.1, 1, 0.1 , AxesLabel
"t min ", "T t
C "
T t
100
C
80
60
40
20
0
5
10
15
20
25
30
t min
36
Ordinära differentialekvationer och Mathematica
HH/ITE/BN
Exempel: En metallbit placeras för avsvalning under rinnande vatten
med temperaturen 10 C. Efter 20 s är temperaturen i metallbiten 150 C
och efter 45 s är den 80 C. Hur varmt är metallstycket två minuter efter
det att man inledde avsvalningen och hur varmt var det då ? Antag att
Newtons avsvalningslag gäller.
Lösningsförslag: Newtons avsvalningslag kan lösas antingen som linjär eller separabel (ODE). Vi tar med första inspektionstillfället för att bestämma integrationskonstanten C1 .
TAvt
DSolve
T t
140
T' t
k t 20
k 10
T t
, T 20
150 , T t , t
Simplify
10
Slutligen fixerar de återstående två inspektionstillfällena proportionalitetskonstanten k samt svaren på frågorna.
TÅk
Solve T t
k
0.0277259, T 0
80., T 0 , T 120
253.754, T 120
. TAvt . t
45, 0, 2
60
18.75
En liten svalkande bild !
PlotT t
T t
. TAvt . TÅk, t, 0, 300 , PlotStyle
Red, AxesLabel
"t
s ", "T t
C " 
C
120
100
80
60
40
20
50
100
150
200
250
300
t s
Exempel: En tank i form av en rak cirkulär kon med spetsen vänd nedåt är helt fylld med
vatten. Så öppnas en kran i botten så att vattnet strömmar ut med en hastighet som i varje
ögonblick är proportionell mot kvadratroten ur vattendjupet, Torricellis lag. Hur lång tid
tar det att tömma tanken om den efter T s är tömd till hälften?
1
 Πr2h
t 3
Lösningsförslag: De två första meningarna ger oss modellen
k h , där h är vattendjupet betraktat som en koordinat
för vattenytans läge i förhållande till konens spets, r är vattenytans radie vid h och k proportionalitetskonstanten. Men om konen har
toppradien R och höjden H ger likformiga trianglar att
1
R
 Π H
t 3
2
h h
r
R
h
,
H
1
R 2
 Π H  h3 
t 3
k h
så
R 2
Π H  h t 2 h ' t
k h
k ht
ht
32
h' t
K.
Vi har nu att beakta
h t 3 2 h' t
h0 H
BVP
hT
H
2
K
ODE
BV
RV
Här är (ODE) separabel. Med BV och (RV) som integrationsgränser möbleras ett ekvationssystem som bestämmer K och tömningstiden Τ
H 2
ekv

T
h3
2
h
H
1


20
2
0
3
 K t,  h
0
8 H 5 2
H
2
H5 2
K T,
K Τ
5
Τ
2
h
 K t
0
HH/ITE/BN
Ordinära differentialekvationer och Mathematica
37
Solve ekv, K, Τ
N
2 H5 2
8 H5 2
K
4
,Τ
2 T

20 T
8 T
31
0.329289 H 5 2
K
,Τ
1.21474 T
T
Eller en riktigt brutal version som går direkt på "urmodellen" och bestämmer såväl h t som "ursprungligt" k.
1
hÅk
2
R
DSolveD
h t
Π
3
h t , t
k t
h t
,
H
H
k' t
0, h 0
H, h T
, h t , k t
, t
2
H 5 2  2 t 8 t 8 T
2
2
H  2 ΠR
8ΠR 
k t
25
5
2
8Π
R2 
20 T
25
H  2 Π
8Π
R2 
2
5
2
20 T
2
25
T
,h t
, k t
20 T
5
H 5 2  2 t 8 t 8 T
T
,h t
,
2
R2
25
T
,h t
2
H 5 2  2 t 8 t 8 T
H  2 Π
8ΠR 
, k t
20 T
k t
2
H  2 ΠR
T
,h t
R2
H 5 2  2 t 8 t 8 T
2

2
5
2
Här duger bara den första lösningen, med tömningstid som ovan.
Solveh t
Τ
t
0 . hÅk 2 , t
N
First
1.21474 T
Till slut vattenytans reseberättelse mot botten
Ploth t
. hÅk 1
PlotRange
.H
16 . T
1, t, 0, t
0, 16 , AxesLabel
"t T
T . Τ , PlotStyle
s ", "h t
Blue,
cm " 
h t cm
15
10
5
0.0
0.2
0.4
0.6
0.8
1.0
1.2
tT s
Exempel: En fisk med massan m slutar simma vid hastigheten v0 och glider
horisontellt genom vattnet. Bestäm fiskens hastighet som funktion av läget
samt hastighet och läge som funktion av tiden om friktionskraften mot vattnet
är proportionell mot hastigheten med proportionalitetskonstanten k. Hur långt
rör den sig ? När stannar den?
Lösningsförslag: Vi känner igen Newton, m x F. Vi börjar med att rulla ut måttbandet, det vill säga x-axeln, precis då den slutar
simma. Den enda kraft som verkar i färdriktningen är då motståndet från vattnet k x. Vi får då
BVP
mx
x0
kx
0
x0
v0
ODE
BV
I första fallet är vi primärt intresserade av hastighet som funktion av läget x x . Vi börjar med den typiska handräkningsversionen
x
t KR
som man så ofta hittar i Mekanikböcker och gör därför omskrivningen x
x
x
x
t
x
x
.
x
Därmed är vi över i en separabel
(ODE) som vi integrerar direkt med (BV) som gränser
mx
x
x
kx
x
v0
x
x
0
k
m
x
x
x
v0

x
k
x
m 0
x
v0
k
x.
m
38
Ordinära differentialekvationer och Mathematica
HH/ITE/BN
Så hastigheten avtar linjärt med läget eller sträckan den glidit. För att få hastighet och läge som funktion av tiden kan vi antingen ge
k
x
m
v0
oss på den nyss erhållna x
eller den ursprungliga m x
k x. Notera att i båda fallen är x
x t . Vi provar båda
vägarna . Den förstnämda är både separabel och linjär. Eftersom vi har haft separationsångest ett tag så väljer vi att öva lite på en
linjär av första ordningen
k
k
x
m
v0
x
x
k
x
m
v0
k
m
Slutligen BV x 0
t
k
mv0
k
x
k
mv0
k
0: 0

IF
C1
m
m
k
t
m
m
t
C1
0
C1
t
k

t

m
mv0
k
x
m
k
mv0
k
x
k
t
C1
m
t
v0
Separabel
t
lösningen till BVP är x t
k
mv0
1
k
m
t

Hastigheten får vi genom att derivera
x
t
x

t
mv0
1
k
k
t
m

mv0
0
k
k
k

m

m
t
k

v0
m
t
Så hastigheten avtar exponentiellt med tiden den glidit. Man ska alltid hålla vad man lovat varför det nu är dax för den andra vägen
med ursprungliga (ODE) som är en linjär andra ordningens homogen med konstanta koefficienter.
kx
mx
k
2m
karakteristiska ekvationen och dess två rötter r1,2
k
x
m
x
k
0
2
2m
r1
0
0
r2
k
k
m
Fall 1
C2
m
k
Eftersom ODE är homogen så är x t
xh t
xp t
C1
xh t
C1
C2
m
t
t
0
k
Slutligen BV 
x0
0
x0
v0
C1
:
C2
k
m
mv0
k
Glidsträckan blir ett gränsvärde x t
m
k
C2
m
då t
0
0
C1
0
C2
v0
mv0
k
mv0
k
återigen lösningen till BVP x t
DSolve
m x '' t
k x' t , x 0
0, x ' 0
v0 , x t , t
First
m


x t
k
PlotEvaluate
PlotStyle
xAvt 1, 2 , D xAvt 1, 2 , t
Red, Blue , PlotRange
xt ,xt
1.0
0.8
0.6
0.4
0.2
1
2
3
4
5
t
Exempel: En ögonblicksbild av en höjdhoppare precis vid
upphoppet kan beskådas i figuren. Sök v0 och Θ så att
hopparens tyngdpunkt G precis klarar höjden.
. m
1, k
All, AxesLabel
1, v0
1
k
v0
Simplify
kt
m v0 1
k
m
t

. Så resan är begränsad i rummet men inte i tiden. Den stackarn glider för gott!
Avslutningsvis kollar vi om det moderna lösningssättet med Mathematica delar vår åsikt. Därefter väljer vi m
kvalitativ överblick av x t och x t under resan.
xAvt
mv0
1
k
, t, 0, 5 ,
"t", "x t , x t "
1 och får en
HH/ITE/BN
Ordinära differentialekvationer och Mathematica
39
Lösningsförslag: Placera ett koordinatsystem med origo där G befinner sig enligt figur. Ta sedan hjälp av Newton och ställ upp
rörelseekvationerna i x- och y-riktningarna. Applicera begynnelsevärden och lös begynnelsevärdesproblemet.
xÅy
DSolve
m x '' t
m y '' t
0, x 0
0, x ' 0
v0 Cos Θ ,
m g, y 0
0, y ' 0
v0 Sin Θ
1
x t
t v0 cos Θ , y t
2
2 t v0 sin Θ
, x t ,y t
,t
First
g t2 
När klockan är t har vi nått högsta punkten på banan och då är y 0. Detta möblerar ett ekvationssystem vars lösning är den önskade
informationen. Vi har att uppfylla tre villkor (ekvationer) med tre obekanta, v0 , Θ och stigtiden t som bonus.
Solve x t
v0
v0
1, y t
5.04228, Θ
5.04228, Θ
1.06, y ' t
1.13005, t
2.01155, t
0
0.464872 , v0
0.464872 , v0
. xÅy . DxÅy, t . g
5.04228, Θ
5.04228, Θ
9.81, v0 , Θ, t 
2.01155, t 0.464872 ,
1.13005, t 0.464872
Eftersom vi inte befattar oss med negativa farter, tider eller vinklar är det bara den sista lösningen som duger.
I höjdhopp behöver inte tyngdpunkten passera ovanför ribban. Detta var svårt att uppnå på "Hedenhös" tid då man använde den så
kallade saxstilen. Lite bättre blev det med dykstilen och dess utveckling med så kallat "hängande ben", men det verkliga genombrottet kom under OS 1968 i Mexico City då amerikanen Dick Fosbury (1947-), på bild nedan, chockade friidrottsvärlden med att visa
upp en helt ny hoppstil, "the Fosbury Flop", som han vann guld med på 2.24. För den riktigt vige erbjuder denna stil just en rejäl
sänkning av tyngdpunkten under det att kroppen krånglar sig över ribban. Vi som känner till mekanikens energilagar vet att detta
betyder ökad hopphöjd för given utgångsfart.
Exempel: Ungefär 25 km uppströms Fylleån, den som rinner strax söder om Eurostop,
ligger i serie fyra sjöar som väsentligen genomsköljs av nämnda å. Dessa är i ordning
Gyltigesjön Töddesjön Simlången Brearedssjön
6
3
Volym 10 m 
3.4
1.2
7
1.2
I figuren till höger ses Töddesjön med sin största ö Storön. Antag att det under en kort
tid sker ett utsläpp av en förorening på en plats uppströms Gyltigesjön. Studera hur
denna vandrar genom sjöarna för att så småningom klinga av. Flödet i Fylleån antar
vi vara 7000 m3 dygn och att sjöarna är helt rena från början.
Lösningsförslag: Antag att utsläppet är
Plott
t
, t, 0, 10 , PlotStyle
Orange, AxesLabel
"t
dygn ", "c0 t
mg m3 "
c0 t mg m3
0.3
0.2
0.1
t dygn
2
4
6
8
10
Vi har typiskt blandningsproblem med perfekt omrörning enligt tidigare beskrivning. Men med fyra obekanta koncentrationer som
möblerar ett system med fyra (ODE) är det helt uteslutet att räkna för hand. Vi tar istället hjälp av vår gamle arbetshäst NDSolve
och löser det kopplade differentialekvationssystemet.
40
Ordinära differentialekvationer och Mathematica
q 7000;
c0 t t ;
Tstud 10 365;
cAvt
HH/ITE/BN
Studera systemet under tio år
6
NDSolve3.4
10 cG ' t
q c0
1.2
106 cT ' t
q cG t
cT t
,
Töddesjön
7.0
106 cS ' t
q cT t
cS t
,
Simlången
cG t
,
Gyltigesjön
1.2 106 cB ' t
q cS t
cB t ,
Brearedssjön
cG 0
cT 0
cS 0
cB 0
0
BV
First;
, cG t , cT t , cS t , cB t , t, 0, Tstud 
Å så här ser skådespelet ut under de tio första åren
PlotEvaluate Last
PlotRange
cG t ,cT t ,cS t ,cB t
cAvt . t
All, AxesLabel
365 t , t, 0, 10 , PlotStyle
"t
Orange, Green, Blue, Red ,
år ", "cG t ,cT t ,cS t ,cB t
mg m3 "
mg m3
0.0020
0.0015
0.0010
0.0005
2
4
6
8
10
t år
Fotnot: Invid Gyltigesjön ligger det som var Halmstadbornas vintersportparadis från 1950-75. Bestyckade med träskidor, råttfällor,
anorak, mössa och lovikavantar begicks mången debut i utförsåkningens ädla konst. Inspirerade av vad man såg i den i hemmen allt
vanligare så kallade TV-apparaten drog man varje snörik helg till "Brearedsfjällen". Med dagens "bortskämda" referenser är de
närmast att betrakta som "vägbulor".
Anstormningen av rejält nedfrusna frisksportare medförde att Turistfrämjandet lät uppföra en sportanläggning, Simlångsgården.
Släplift på egen risk installerades och längdskidspår anlades. Nu blev det riktigt folkligt. För de som inte tyckte sig ha råd att köpa
fika fanns värmestugor där medhavda mackor och varm chocklad (med skinn!) kunde inmundigas under livliga diskussioner kring
teknikgrenarna. Med charterresornas uppgång kom dock verksamheten i svårigheter. Men enligt senaste rapporter är nu verksamheten på frammarsch och i backarna har åtminstone pulkafolket hittat tillbaka medan de gamle myser och ventilerar minnen från
svunna tider, kan ett samtida ögonvittne från både då och nu berätta.
Exempel: Volterra–Lotkas jägar–byte modell. Studera
ett stort skogsparti där endast rävar och kaniner lever.
Uppgiften blir att studera djurens utveckling över tiden. När de lever i den bästa av världar, det villl säga ostörda, bygger modellen
på följande förutsättningar
Vid tidpunkten t finns det k t st kaniner och r t st rävar.
Rävarnas enda föda är kaniner vars föda det finns i överflöd.
Vid avsaknad av rävar gör kaniner det dom är bäst på och ökar sin population enligt Malthus lag, k ' t
Vid avsaknad av kaniner svälter rävarna ihjäl. Även detta enligt Malthus lag, r' t
Αk t .
Βr t .
Då det finns både kaniner och rävar i skogen är det rimligt att anta att antalet kaniner som dödas per tidsenhet av rävarna är
proportionellt mot både antalet kaniner och antalet rävar. Kaninekvationen korrigeras alltså till k ' t Αk t Γ k t r t .
Med samma betraktelsesätt får vi då att rävarnas antal ökar. Så rävekvationen korrigeras till r' t
Då har vi
Βr t
Δk t r t .
HH/ITE/BN
Ordinära differentialekvationer och Mathematica
k' t
r' t
k0
r0
BVP
Αk t Γ k t r t
Βr t Δk t r t
k0
r0
41
ODE
BV
Detta är ett system av icke-linjära differentialekvationer av första ordningen. De går inte att lösa analytiskt utan vi får tillgripa den
numeriska lösaren i Mathematica. Vi väljer lite numeriska värden och öppnar skådespelet.
kÅr t
k t ,r t
. NDSolve
k' t
4k t
0.02 k t r t ,
r' t
0.9 r t
0.001 k t r t ,
k 0
200, r 0
200 , k t , r t
t, 0, 25
First;
,
Några bilder kan inte skada
PlotEvaluatekÅr t , t, 0, 25 , PlotStyle
Gray, Red , AxesLabel
"t", "kaniner,rävar" 
kaniner,rävar
2500
2000
1500
1000
500
5
10
15
20
25
t
Man kan rita lösningarna på parameterform och bilda ett så kallat fasporträtt. Då t ökar vandrar punkten k t , r t med t som
parameter utmed en sluten kurva enligt pilen i figuren nedan. Om vi t.ex. befinner oss i punkten A så kan vi förvänta oss att i en nära
framtid så ökar antalet rävar medan antalet kaniner minskar.
ParametricPlotEvaluatekÅr t , t, 0, 5 , PlotRange
AxesLabel
"kaniner", "rävar" , AspectRatio
0, 2500 , 100, 350
0.6, Epilog
, PlotStyle
Red,
Red, Arrowheads 0.06 ,
ArrowkÅr 2.2 , kÅr 2.21 , GrayLevel 0 , Text"A", kÅr 2.1 , Background
White
rävar
350
300
250
A
200
150
100
0
500
1000
1500
kaniner
2500
2000
Vi kan göra en ekologisk betraktelse genom att studera en bukett av kurvor sprungna ur olika (BV). Antag att vi befinner oss i punkt
P1 i figuren nedan och att vi efter t.ex. ett mänskligt ingripande "reducerar" antalet kaniner så vi hamnar i punkten P2 på en annan
kurva i fasplanet. En tid senare får detta som effekt att vi befinner oss i Q2 i stället för Q1 , det vill säga antalet kaniner är många flera
än vad de skulle varit om vi lämnat skogen ifred. Bra exempel på kontraintuitiv mänsklig insats!
rävar
350
300
250
200
P2
Q1
P1
Q2
150
100
0
500
1000
1500
2000
kaniner
2500
Denna modell introducerades av den italienske matematiken Vito Volterra i samband med hans studier av fiskbeståndet i Adriatiska
havet under 1920-talet. Modellen är generell och kan genom olika val av proportionalitetskonstanter och (BV) anpassas till olika
arter av djur.
42
Ordinära differentialekvationer och Mathematica
HH/ITE/BN
Exempel: Svängningar. Vår omvärld är full av företeelser som karakteriseras av svängningar, exempelvis vibrationer i en maskin,
stötar i ett knä vid löpning eller brus i en mobiltelefon. Många av dessa kan med fördel lineariseras till en diskret modell med en
massa, en fjäder och en dämpare. Inte sällan ser man flera sådana "atomer" sammankopplade i större strukturer. Med hjälp av
Newtons rörelselagar får vi så en modell som bör finnas i verktygslådan hos varje praktiserande ingenjör.
För fjädern gäller att kraften som krävs för att förlänga den en sträcka x är proportionell mot denna, det vill säga F f
kx, där k
kallas fjäderkonstanten med enheten N/m. Observera att uttrycket gäller med tecken på x, det vill säga även då den trycks ihop.
Ff
k
Ff
x
x
För dämparen gäller däremot att kraften är proportionell mot förlängningshastigheten (farten), det vill säga Fd cx , där c kallas
dämpningskonstanten med enheten Ns/m. En dämpare består i princip av en cylinder fylld med olja. En rörlig bricka med hål i
delar in cylindern i två kammare. Oljans viskositet, det vill säga hur trögflytande den är, antal hål i brickan och deras storlek avgör
hur trögt det är att flytta den, eftersom olja från ena kammaren då skall flyttas till den andra. Detta ger dämparen dess karakteristiska
funktion som också brukar kallas viskös dämpning. Kraften som krävs för att flytta på brickan är alltså inte beroende av dess läge
utan endast på hastigheten (farten). Jämför potatispress!
Fd
c
Fd
x
x
Massa och fjäder utan yttre last (fri odämpad svängning)
Vi börjar med en massa kopplad till en fjäder. Vi frilägger och
noterar att den motriktade kraften i fjädern är proportionell mot
läget med fjäderkonstanten k. Formulera nu rörelseekvationen
med hjälp av Newton
mx
kx
mx
kx
0
k
x
m
x
0
k
m
Vi känner igen en linjär andra ordningens (ODE) med konstanta koefficienter. Här kallas Ωe
enheten [rad/s]. Man pratar också om egenfrekvensen fe
Ωe
2Π
och egensvängningstiden Te
1
.
fe
för egenvinkelfrekvens med
Prefixet "egen" kommer av att Ωe
endast beror på systemets "egna" inneboende parametrar. Lösningen är en så kallad harmonisk svängning (sinusformad) som fortgår
i all evighet om systemet störs från sitt jämviktsläge. Notera speciellt var Ωe dyker upp!
DSolve m x '' t
kx t , x t , t
k t
x t
c2 sin
k t
c1 cos
m

m
Som väntat får vi två konstanter som fixeras av (BV). Uttrycket ovan brukar ofta skrivas om som en ren sinus- eller cosinusvåg.
a sin Ωt
a2
b cos Ωt
a
b2
b
sin Ωt
a2 b2
och identifierar amplitud R
och fasvinkeln
a2
a2
b2 cos
a2
b2 sin Ωt
sin Ωt
b2 och fasvinkel
cos Ωt
a2 b2
sin
cos Ωt
1
a
b
1, 1
a2 b2
1
a2 b2
Summaformel för sinus
a
som naturligtvis mäts i radianer. Punkten  R ,
bestäms likt argumentet för ett komplext tal. Var och en av de oändligt många vinklar
b

R
ligger på enhetscirkeln
som löser ekvationerna
HH/ITE/BN
Ordinära differentialekvationer och Mathematica
a
b
43
R cos
R sin
kan göra anspråk på att kallas för fasvinkel. På grund av periodiciteten hos cosinus och sinus skiljer de sig åt med en multipel av 2Π
så alla ger "samma avtryck" i enhetscirkeln. Ofta nöjer man sig med den så kallade principalvinkeln som ligger i intervallet Π, Π .
När man räknar för hand gäller det att se till så man hamnar i rätt kvadrant!! Om a 0 eller b 0 är det ju enkelt annars går det bra
b
arctan a 
att beräkna fasvinkeln som
Π eftersom arctan() levererar vinklar i första och fjärde kvadranten. Den avslutande
om a 0
korrektionen kommer sig naturligtvis av att vi kan ha dividerat bort "negativ" information, ty
b
a
b
a
och
b
a
b
.
a
I Mathematica är
det inga problem, speciellt inte om vi använder versionen med två argument ArcTan[a,b]som alltid levererar rätt vinkel i
intervallet Π, Π och givetvis i radianer. Vi tar ett exempelvis med m 2 kg, k 18 N m och BV x 0 1, x 0 1.
xAvt
DSolve
2 x '' t
18 x t , x 0
1, x ' 0
1 ,x t ,t
First
1
x t
sin 3 t
3 cos 3 t 
3
Plot x t
. xAvt, t, 0, 10 , PlotStyle
Red, AxesLabel
"t", "x t "
xt
1.0
0.5
2
4
6
8
10
t
0.5
1.0
a, b
Coefficient x t
. xAvt, Sin 3 t , Cos 3 t
1
 , 1
3
a2
R
b2
10
3
ArcTan a, b
tan
1
3
En sista ängslig koll...
R Sin 3 t
x t
. xAvt
Simplify
0
Massa, fjäder och dämpare utan yttre last (fri odämpad svängning)
Massa kopplad till en fjäder och dämpare. Vi frilägger och noterar att den
motriktade kraften i fjädern är proportionell mot läget och i dämparen mot
hastigheten med dämpningskonstanten c. Formulera rörelseekvationen
med hjälp av Newton
mx
kx
cx
mx
kx
cx
0
x
c
x
m
k
x
m
0
Vi känner igen en linjär andra ordningens (ODE) med konstanta koefficienter. Beroende på om rötterna till den karakteristiska
ekvationen är reella och olika, reella och lika respektive komplexkonjugerade använder vi beteckningen överkritisk, kritiskt
respektive underkritisk dämpning. Detta avgörs av den nya aktören, dämparen, som i praktiken bidrar med att "äta energi", vilket
ger sig tillkänna som en exponentiellt avtagande amplitud hos lösningen. Även den tidigare nämnda egenvinkelfrekvensen Ωe
c
kommer med i lösningen på "samma ställe" som tidigare. Man brukar också definiera dämpningsfaktorn Ζ genom 2Ωe Ζ m .
44
Ordinära differentialekvationer och Mathematica
För de tre fallen ovan gäller Ζ
följande figurer.
1, Ζ
0 respektive Ζ
1. De tre möjliga fallen av homogena lösningar illustreras kvalitativt i
xt
xt
HH/ITE/BN
xt
t
t
t
Fall 2: Kritisk dämpning.
Lösningskurvan skär x axeln högst
en gång. Typiskt utseende för en ny
hjulupphängning på en bil.
Fall 1: Överkritisk dämpning.
Inte så vanligt förekommande i
ingenjörssammanhang.
Exempelvis underkritisk dämpning med m
xAvt
DSolve
1
x t
t6
3 x '' t
3 kg, c
5x t
1, x ' 0
1 ,x t ,t
59 cos
6
ab.ab
. ab
First

Coefficientx t
59 t
, Cos
. xAvt, Sin
6
3
1.
6
59 t
R
1 och x 0
59 t
59 sin
59
5 N m med (BV) x 0
1 x' t , x 0
59 t
7
1 Ns m, k
Fall 3: Underkritisk dämpning.
Mycket vanligt insvängningsförlopp.
Typiskt utseende för en gammal
hjulupphängning på en bil.

6
t3
6
59
Plot R, x t
. xAvt, R , t, 0, 20 , PlotStyle
PlotRange All, AxesLabel
"t", "x t , R t "
Blue, Red ,
xt , Rt
1.5
1.0
0.5
5
10
15
20
t
0.5
1.0
1.5
Massa, fjäder och dämpare med yttre last (påtvingad svängning)
Avslutningsvis kan man naturligtvis utsätta de båda modellerna ovan för en yttre,
ofta kallad störande eller påtvingad, kraft. I nästan alla intressanta applikationer
är denna av periodisk natur F t F0 sin Ωt . Vi väljer att exemplifiera med en
full massa–fjäder–dämpare–modell. Formulera nu rörelseekvationen med
hjälp av Newton
mx
kx
cx
Ft
mx
kx
cx
Ft
x
c
x
m
k
x
m
Ft
m
Vi känner till att lösningen delas upp i "systemets inneboende" homogena lösning xh samt en partikulärlösning x p beroende på
högerledet, ofta kallad fortvarighetslösning (eng. steady-state solution). Allmänna lösningen blir sedan x
xh
x p . Homogena
lösningen brukar kallas för transienta lösningen (eng. transient solution) eftersom den klingar av ganska snabbt med tiden. Kvar
blir x p , därav namnet. När vi löser en differentialekvation med Mathematica får vi inte denna, ibland önskade, uppdelningen. Ett
trick är att lösa systemet utan begynnelsevärden. De termer som då har "konstanterna" Ci på sig tillhör den homogena lösningen och
öppnar därmed för en separation av lösningen.
HH/ITE/BN
Ordinära differentialekvationer och Mathematica
Exempelvis med m
(BVP).
xBVP
1 kg, c
DSolve
1
t
x t
20
2 Ns m, k
5 N m, F0
x '' t
2 x' t
x 0
0, x ' 0
t
37 sin 2 t
80 
10 N, Ω
2 rad s med (BV) x 0
5x t
10 Sin 2 t ,
1 ,x t ,t
First
t
45
0 och x 0
1. Först hela lösningen till
Simplify
1 cos 2 t 
34
Plot x t
. xBVP, t, 0, 15 , PlotStyle
Red, AxesLabel
"t", "x t "
xt
2
1
2
4
6
8
10
12
t
14
1
2
Nu över till separationen x
xODE
1
t
17
0, c2
xp
x p . Samma som ovan men utan begynnelsevärden (BV).
DSolve x '' t
x t
Så c1
xh
2 x' t
10 t  sin 2 t
17 c1
5x t
17 c2
10 Sin 2 t , x t , t
First
Simplify
40 t  cos 2 t 
0 ger oss partikulärlösningen x p och sedan xh
x t
. xODE . C 1
0, C 2
0
Simplify
10
4 cos 2 t
sin 2 t
17
xh
x t
1
t
. xBVP
37 sin 2 t
xp
Simplify
80 cos 2 t
34
Nu kan vi inspektera de olika delarna. Allmänna lösningen anpassar sig raskt till partikulärlösningen eftersom den homogena
lösningen snabbt klingar av med tiden beroende på exponentialfunktioner med negativa exponenter.
Plot x t
. xBVP, xh, xp , t, 0, 15 ,
PlotStyle
Red, Blue, Green , AxesLabel
"t", "x t ,xh t ,xp t "
x t ,xh t ,x p t
2
1
2
4
6
8
10
12
14
t
1
2
Även här brukar fortvarighetslösningen skrivas om till ren sinusvåg x p t
vinkelfrekvens Ω, här = 2, och fasvinkel
a, b
10

Coefficient xp, Sin 2 t , Cos 2 t
40
,
17

17
a2
R
10
17
i radianer naturligtvis.
b2
asin Ωt
bcos Ωt
Rsin Ωt
med amplitud R,
46
Ordinära differentialekvationer och Mathematica
HH/ITE/BN
ArcTan a, b
tan
1
4
En sista ängslig koll...
R Sin 2 t
xp
Simplify
0
Amplitud över alla gränser
En enkel massa med fjäder och dämpare är uppriggat
enligt figur. Sök fortvarighetslösningens amplitud
orsakad av den störande kraften.
k
m
Om den störande kraftens vinkelfrekvens Ω närmar sig egenvinkelfrekvensen Ωe
kommer lösningens amplitud att öka
dramatiskt. Man talar om resonans (eng. resonance). Vid konstruktion med liten dämpning ska man alltså undvika detta tillstånd
eftersom ett haveri är att vänta. Vi illustrerar med en bukett dämpare samt Ωe 1 och ritar upp partikulärlösningens amplitud R.
amp
;x t
Dor
: A Cos Ω t
B Sin Ω t
Solve Coefficient x '' t
AppendTo amp, Norm
3
5
, 1,
;
, c,
10
100
A, B
c x' t
x t
. First r
Sin Ω t , Cos Ω t , Sin Ω t
0, A, B
;
;
Plot Evaluate amp , Ω, 0, 2 , PlotStyle Hue
Range 1, 1
PlotRange
0, 2 , 0, 4 , AxesLabel
"Ω Ωe ", "R"
15,
1
15 ,
R
4
3
2
1
0
0.0
0.5
1.0
1.5
2.0
Ω Ωe
Vi poängterar ytterligare de två fundamentala och viktiga fallen då den störande vinkelfrekvensen Ω
extra tydlig om vi illustrerar utan dämpare.
Fallet Ω
Ωe för modell utan dämpare
Speciella tillämpningar har vi då Ωe
xAvt
DSolve
Ω är mycket mindre än Ωe
x '' t
36 x t
x 0
0, x ' 0
Sin 5 t ,
1 ,x t ,t
sin 5 t
sin 6 t 
11
Denna kan, med ett litet trick, skrivas om till
TrigFactor x t
2
11 t
11
t
cos
sin
2
2
. xAvt . 6
.
Ω, så exempelvis med Ωe
First
1
x t
xt
Ωe och Ω
6
Simplify
6 och Ω
5.
Ωe . Effekten blir
HH/ITE/BN
Ordinära differentialekvationer och Mathematica
2
t
cos 2 
11
Här kan vi nu betrakta
47
som en med tiden "sakta" varierande amplitud till den "fortare" svängande sinusvågen sin
11 t
.
2
Bland musiker brukar detta kallas för svävning (eng. beat) och kan höras strax innan en gitarrsträng är riktigt stämd. Elektroingenjörer identifierar fenomenet som amplitudmodulering (eng. amplitude modulation) och hittas på radions AM-band. Där känner vi
också igen FM-bandet baserat på frekvensmodulering (eng. frequency modulation).
2
t
Cos ,
11
2
Plotxt,
t
Cos , t, 0, 6 Π , PlotStyle
11
2
2
t
"t", "x t ,
cos
"
11
2
AxesLabel
2
2
Red, Blue, Blue ,
t
cos
xt ,
11
0.2
2
0.1
5
10
t
15
0.1
0.2
Fallet Ω
Ωe för modell utan dämpare
Detta resonansfall måste ovillkorligen undvikas vid all konstruktion eftersom det leder till haveri. Samma som föregående med
Ωe Ω 6.
xAvt
DSolve
13
x t
x '' t
36 x t
x 0
0, x ' 0
Sin 6 t ,
1 ,x t ,t
First
Simplify
Apart
1
sin 6 t
72
t cos 6 t 
12
Här kommer nu partikulärlösningens tidsberoende amplitud
allt mera kraftiga svängningar eftersom
t
Plotx t
. xAvt,
t
12
att snabbt dominera förloppet. Om inget hejdar utvecklingen får vi
, med ett oundvikligt haveri som följd.
t
,
12
då t
t
12
, t, 0, 6 Π , PlotStyle
Red, Blue, Blue ,
12
t
AxesLabel
"
"t", "x t ,
12
t
xt ,
12
1.5
1.0
0.5
5
10
15
t
0.5
1.0
1.5
Ett gammalt välkänt exempel på denna typ av haveri är kollapsen av Tacoma Narrows Bridge 1940, där den måttliga vindens
frekvens matchade längden på bron och satte en torsionsmod i resonanssvängning. På www hittar man lätt filmsnuttar från skådespelet! Ett exempel från senare tid är rymdfärjan Challenger där en turbopump oavsiktligt konstruerades med resonans. Lyckligtvis
upptäcktes det före start, men felet kostade NASA flera miljoner dollar. Vid olyckan 1985 kom problemställningen att diskuteras på
nytt.
48
Ordinära differentialekvationer och Mathematica
HH/ITE/BN
Exempel: Ett så kallat bungy jump är väl bekant för de flesta. Från
en hoppställning högt ovan mark kastar man sig handlöst ut. Livlinan
utgörs av en gummisnodd vars ena ände är fäst vid uthoppsplatsen
och den andra surrad runt vristerna på offret. Välj lite numeriska
data och gör en simulering av ett bungy jump
Lösningsförslag: Om vi placerar en y-axel med origo vid marken och pekande uppåt har vi återigen Newton, my F. Krafterna
som verkar är tyngdkraften och den från gummisnodden. Eftersom en gummisnodd, för enkelhets skull antar vi, fungerar ungefär
som en fjäder frånsett att den inte kan ta trycklaster, är den lite kinkig att simulera och lösa analytiskt. Men med storsläggan
NDSolve går det både snabbt och enkelt. För att offret inte ska svänga upp och ned i all oändlighet, vilket är fallet med en fjäder,
lägger vi in lite dämpning i snodden. Med fjäderkonstanten k, dämpning c och den naturliga längden L gäller om hoppställningen har
höjden H att kraften i snodden
Fgs
kH
0
y
L
c y om H
om H
y
y
L
L
0
0
Vi får då
BVP
my
mg
y0
H
y0
0
Fgs
ODE
BV
Nu gäller det bara att översätta den styckvis definierade funktionen Fgs till Mathematica, vilket lätt görs med funktionen Piecewise eller If. Vidare kräver NDSolve att resan görs helt numeriskt, så välj t.ex. m
c 50 Ns m, så har vi direkt en studie över de T första sekunderna
T
75 kg, H
15;
Fgs

yAvt
50 m, L
30 m, k
125 N m och
Restid
k H
y t
NDSolve
L
0
c y' t
H
H
y t
y t
L
L
0
;
0
F i gummisnodden
m y '' t
m g Fgs ,
y 0
H, y ' 0
0
. m 75, g 9.81, H 50, L
y t , t, 0, T
First;
30, k
125, c
50 ,
ODE
BV
Data
En bild över spektaklet piggar alltid upp
PlotEvaluate
PlotStyle
y t
. yAvt, D y t
. yAvt, t
Red, Blue , AxesLabel
"t
, t, 0, T ,
s ", "y t
m , y t
m s "
yt m,yt ms
40
20
2
4
6
8
10
12
14
t s
20
Man ser tydligt på hastighetsgrafen när gummisnodden börjar bromsa in det fria fallet.
Exempel: En fallskärmshoppare, vars vikt är 75 kg, hoppar ut
från en helikopter 4 km ovan markytan. Antag att luftmotståndet
är proportionellt mot hopparens fart. Proportionalitetsfaktorn är
15 kg s när fallskärmen är outlöst och 105 kg s när fallskärmen
är utvecklad. Antag att fallskärmen utvecklas 1 min efter uthoppet
från helikoptern. Sök läge och hastighet under resan samt bestäm
hur lång tid hoppet tar.
Lösningsförslag: Om vi placerar en y-axel med origo vid marken och pekande uppåt har vi återigen Newton, m y F. Krafterna
som verkar är tyngdkraften och den från fallskärmen. Vi har att ta ställning till två stycken (BVP). Ett under den första minuten med
HH/ITE/BN
Ordinära differentialekvationer och Mathematica
kända (BV), sedan ett annat under resten av nedfärden. Denna måste matchas med BV
att den sammanfogas med lösningen under den första minuten. Eller lite mera precist
BVP1
m y1
y1 0
mg k1 y1
H
y1 0
0
ODE
BV
för t
0, 60
och
BVP2
49
RV vad gäller tid, läge och hastighet så
m y2
y2 60
mg k2 y2
y1 60
y2 60
y1 60
ODE
BV
för t
60,
När man så äntligen trött kommer i mål gäller det att mata Plot med rätt lösning i rätt intervall ;-(. Räddaren heter naturligtvis
NDSolve. Vi löser (BVP1) för hela resan t 0,
och byter helt enkelt proportionalitetskonstant k1 k2 vid 60 s on the fly !
Smidigt! Nu är det bara att sätta igång! Kör ordentligt i backen så tar vi reda på efteråt hur lång tid hoppet tog!
yAvt
NDSolve
m y '' t
m g If t Τ, k1 , k2 y ' t ,
y 0
H, y ' 0
0
. m 75, g 9.81, H 4000, Τ 60, k1
y t , t, 0, 500
First;
15, k2
105 ,
ODE
BV
Data
Hopptid
Thopp
t . FindRoot y t
0 . yAvt, t, 200
241.56
Några bilder över spektaklet piggar alltid upp
PrintPlot y t
. yAvt, t, 0, Thopp , PlotRange
All, PlotStyle
AxesLabel
"t s ", "y t
m " ,"
",
PlotEvaluate D y t
. yAvt, t , t, 0, Thopp , PlotRange
PlotStyle
Blue, AxesLabel
"t
s ", "y t
Red,
All,
m s "
yt ms
yt m
4000
50
100
150
200
t s
10
3000
20
2000
30
1000
40
50
100
150
200
t s
50
Man ser tydligt på grafen över hastigheten när fallskärmen börjar bromsa in det fria fallet. Man kan också lätt identifiera de två
gränshastigheterna, det vill säga då vi inte har någon acceleration längre. Dessa blir utan respektive med fallskärm
D y t
. yAvt, t
.t
55, 200
49.0492, 7.00714
Exempel: Hållfasthetslära. I hållfasthetslära och byggstatik böjs det
balkar i tid och otid. Det är ett centralt moment att behärska för ingenjörer
som skall arbeta med konstruktion av exempelvis hus, flygplan, bilar eller
mobilmaster. Därför ingår en kurs i hållfasthetslära på i princip alla
ingenjörsutbildningar. Alla viktiga aspekter såsom utböjning och om
materialet i balken håller för den pålagda lasten går tillbaka på att
bestämma den så kallade elastiska linjen.
Lösningsförslag: Om utböjningen i balken är u x , E x materialets elasticitetsmodul, I x balktvärsnittets yttröghetsmoment och
q x den utbredda lasten, riktad som u x , gäller det att lösa en fjärde ordningens differentialekvation.
E x I x u '' x ''
qx
Denna modell brukar kallas Euler-Bernoulli-balk efter upphovsmännen (1700-talet!) och lösningen u x benämns elastiska linjen.
Generationer av ingenjörsstudenter har "plågats" med ändlöst handarbete för lösa den med två (RV) i varje ände (u utböjning, u '
vinkeländring, u '' moment och u ''' tvärkraft) och några olika vanliga q x . Oftast är EI konstant, så kallad jämnstyv balk. För
Πx
Mathematica är det inga problem och vi exemplifierar med en så kallad konsolbalk enligt figur med q x
q0 sin L .
50
Ordinära differentialekvationer och Mathematica
HH/ITE/BN
Πx
u x
u x
. DSolveEI u '''' x
q0 Sin
,
L
u 0
q0 6 L4 sin

Πx
L
6 Π L3 x

3 Π3 L2 x2
0, u ' 0
0, u '' L
0, u ''' L
0, u x , x
Π3 L x 3 

6 Π4 EI
Först utböjningen sedan moment och tvärsnittskraft.
EI u Ξ L
Plot
, Ξ, 0, 1 , PlotStyle
Orange, AxesLabel
q0 L4
"x L", "q0 L4 EI"
q0 L4 EI
0.2
0.4
0.6
0.8
1.0
x L
0.02
0.04
0.06
EI u '' Ξ L
Plot
EI u ''' Ξ L
,
q0
, Ξ, 0, 1 , PlotStyle
L2
AxesLabel
Blue, Red ,
q0 L
"x L", "q0 L2 EI, q0 L EI"
q0 L2 EI , q0 L EI
0.2
0.2
0.4
0.6
0.8
1.0
x L
0.2
0.4
0.6
Det helt dominerande verktyget idag inom beräkningsteknik för exempelvis hållfasthetslära och strömningsteknik kallas finita
elementmetoden (FEM) och är helt datorbaserad. Om denna kan du läsa i en annan kurs...