NUMERIEKE INTEGRATIE VAN DIFFERENTIAALVERGELIJKINGEN

B.J.G.Hamer augustus 1993

De bedoeling van dit artikel is tweeledig. Sommige differentiaalvergelijken zijn dermate gecompliceerd dat het uiterst lastig, zo niet onmogelijk is deze analytisch op te lossen. Daarnaast komt het voor dat men met differentiaalvergelijkingen geconfronteerd wordt terwijl men daar niet of slecht mee bekend is. Voor bovengenoemde gevallen is dit artikel bedoeld.
In het verleden heb ik de hier onder beschreven methode enkele malen met succes toegepast, het ging in dit geval om differentiaalvergelijkingen die analytisch, althans door mij, niet op te lossen waren. De onderstaande rekenmethode is alleen van toepassing op gewone differentiaalvergelijkingen, geen partieële, en de hoogste afgeleide van de 2e graad. Aan de hand van een simpel voorbeeld zullen we de methode beschrijven.
Iedereen is bekend met de bewegingsvergelijking van Newton:

Dit geldt alleen als de massa niet in de tijd verandert, maar we zullen ons bij aardse snelheden houden.

Voor de versnelling a kunnen we ook schrijven:

Beschouwen we een vrij vallend lichaam dan is a = g en krijgen we dus

 

of

Hier staat dat de linkerterm ongeveer gelijk is aan de rechterterm. de afwijking wordt geringer naarmate de tijdstap,Δt kleiner genomen wordt. De consequentie hiervan is dat de rekentijd toeneemt. De index i betekent dat de hele term waar deze bij staat op het aanvangstijdstip bepaald dient te worden en de term met index i+1 op een tijdstip,Δt later.

We hadden al eerder gezien dat d(dy/dt)/dt = g

zodat

Bij g hoeven we geen index te gebruiken daar deze constant is (voor korte tijd in ieder geval!)

De volgende stap levert een van de af te leiden basisvergelijkingen op.

De term dy/dt is per definitie

De indices hebben hier dezelfde betekenis als bij de boven gegeven vergelijkingen.

Hieruit volgt

Deze twee afgeleide vergelijkingen vormen de basis voor de numerieke oplossingsmethode. De berekening gaat als volgt, met vgl.(2) berekent men uit y0 en dy/dt0 de waarde voor y1. Vervolgens uit dy/dt0 en g de waarde voor dy/dt1. Daarmee is de eerste rekenlus doorlopen en kunnen we deze stappen herhalen met de nieuwe gevonden waarden als uitgangspunt. We stoppen het proces zodra de gewenste eindtijd bereikt is. Het handigst is om deze berekeningen op een PC uit te voeren, maar ook op een programmeerbare rekenmachine gaat het prima, is alleen vaak wat lastiger te programmeren omdat je al gauw in conflict komt met het aantal programmastappen.

Voor het bovenstaande voorbeeld is gekozen omdat dit ook analytisch eenvoudig op te lossen is. Bekijkt men de antwoorden in beide gevallen dan blijkt er een afwijking te zijn. Dit vindt zijn oorzaak in het feit dat we er in het begin van uitgegaan zijn dat de snelheid in het eerste tijdsinterval nul was. We veranderen daarom van strategie, we gaan nu eerst met een halve tijdstap de snelheid berekenen. d.w.z. dy/dt1/2 uit dy/dt0 en g. Vervolgens wordt uit y0 en dy/dt1/2 een waarde voor y1 berekend met een hele tijdstap. De overige berekeningen worden voor zowel y als dy/dt weer met een hele tijdstap gedaan en we krijgen enigszins gewijzigde basisvergelijkingen.

Deze worden

en

Vergelijkt men een berekening verricht met deze 3 vergelijkingen dan blijkt de overeenkomst met de analytische oplossing zeer goed te zijn.

Tot nu toe zijn we uitgegaan van een constante, n.l. g. Het blijkt echter dat de methode even zozeer bruikbaar is voor ingewikkelder vergelijkingen. Daartoe behoeven we alleen de vergelijking in een vorm te schrijven dat hij voldoet aan:

 

Hierbij is C een willekeurige functie van alle grootheden, zoals constanten, y, t, x, dy/dt, dy/dx etc. De vergelijkingen (3) t/m (5) krijgen daardoor de vorm

en

Het is vaak handig om een rekenschema te tekenen, een voorbeeld is gegeven in Fig.1 Dit schema is gemaakt om een volgend voorbeeld toe te lichten uit de mechanica en wel de sterkteleer. Gevraagd wordt de doorbuiging en de doorbuigingshoek te berekenen van een bladveer met rechthoekige doorsnede. Een probleem dat triviaal lijkt omdat dit te vinden is in vrijwel elk handboek. Alleen wordt daarbij stilzwijgend verondersteld dat het gaat om kleine doorbuigingen, meestal niet expliciet aangegeven zodat men behoorlijk de mist in kan gaan. We zullen dan ook zowel het vereenvoudigde geval als het nauwkeurige geval behandelen. Terwille van de overzichtelijkheid gaan we op een andere notatie over.

Een bladveer van elastisch materiaal welke aan de wet van Hooke voldoet zal dusdanig deformeren, hetzij door zijn eigen gewicht, hetzij door een belasting, dat de kromtestraal van de elastische kromme evenredig is met het buigende moment M.