Løsning af differentiallignings-systemer

MATLAB tilbyder flere hjælpefunktioner til numerisk løsning af et system af ordinære differentialligninger, bl.a. ode23 og ode45, der benytter en 2. og 3. hhv 4. og 5. ordens Runge-Kutta-Fehlberg integrations algoritme med (automatisk genererede) variabel skridtlængde. Nærværende afsnit vil koncentrere sig om disse to metoder. Benyt hjælpefunktionerne til at undersøge, hvilde andre integrations algoritmer MATLAB kan tilbyde.

Et system af (ikke-lineære) differentialligninger kan altid omformes til et system af første ordens (ikke-lineære) differentialligninger:

$$\dot{x}=f(x,t)$$

hvor $x = x(t) \in {\bf R}^{n}$ er tilstandsvariablen, $t$ er (sædvanligvis) tiden, og $f$ er en funktion der giver den tidsafledede $\dot{x}$ af $x$ som funktion af $x$ og $t$. Det aktuelle system af (ikke-lineære) første ordens differentialligninger specificeres i en M-fil (se kapitel 3 om M-filer og evt nedestående eksempel).

Med kommandoen

         [t,x] = ode23('name', tspan, xstart);

startes integrationen, af det differentialligningssystem der er specificeret i M-filen ved navn name, ved hjælp af ovennævte 2. og 3. ordens integrations algoritme. Hvis man ønsker at benytte en 4. og 5. ordens integrations algoritme skrives ode45 istedet for ode23. Begge integrations algoritmer skal (mindst) have tre input argumenter. Disse er navnet på en M-fil (som allerede beskrevet), det tidsinterval [tstart tslut], hvori systemet skal løses, og en startbetingelse xstart. Det skal nævnes, at der er mulighed for at specificere yderligere input argumenter, såsom regnenøjagtighed (se MATLABs egen manual for flere detaljer).

Eksempel

Betragt Van der Pol ligningen, der er en anden ordens differentialligning,

$$\ddot{x}+(x^2-1)\dot{x}+x=0$$

Sættes $x_2=x$ og $x_1=\dot{x}$ kan Van der Pol ligningen skrives som følgende system af første ordens differentialligninger:

$$\dot{x}_1 = x_1(1-x_2^2)-x_2$$

$$\dot{x}_2 = x_1$$

Det første skridt i simuleringen er at skabe en M-fil indenholdende disse diferentialligninger. Denne kaldes vdpol.m Dette gøres ved at klikke på File menuen, derefter vælge undermenuen New, og endelig undermenupunktet M-file. I denne M-fil skrives:

     function xdot =  vdpol(t,x)
     xdot          =  zeros(2,1);
     xdot(1)       =  x(1) .* (1 - x(2).^2) - x(2);
     xdot(2)       =  x(1);

I første linie erklæres funktionen og samtidig fastlægges, at det er argumentet til vdpol der returneres. I anden linie nulstilles vektoren xdot samtidig med at den oprettes som en $2\times 1$ matrix. I tredie linie skrives den første differentialligning og i fjerde linie skrives den anden differentialligning. Når indholdet i M-filen er indskrevet i denne skal filen gemmes. Dette sker ved at klikke på File menuen og derefter klikke på undermenuen Gem som ..., herefter åbnes et vindue, hvori filenavnet skrives. I dette eksempel skrives vdpol.m. Da er filen gemt og vinduet kan lukkes. Dette er proceduren for at gemme filer i default editoren 'Notepad'. Ved valg af anden editor, kan proceduren være en anden.

For at simulere Van der pol differentialligningen, som er beskrevet i M-filen vdpol.m over et tidisinterval $0\leq t \leq 20$ vha. ode23 skrives

  tspan     = [0 20];     %Starttid og sluttid 
                          %for integrationen
  xstart    = [0;0.25];   %startbetingelse
  [t,x]     = ode23('vdpol',tstart,tslut,xstart); 
                          %simuleringen udføres kurverne
  plot(t,x)               %x(1) og x(2) plottes versus t 

På figur 4.2 ses resultatet af dette plot.

Figure 4.2: Plot af Van der Pol-systemmet. x(1) og x(2) som funktion af tiden.

De tre første linier giver sig selv. I fjerde linie kaldes ode23 algoritmen, som mindst skal have fire indput argumenter. Bemærk at navnet på den M-fil der indenholder differentialligningerne optræder uden `efternavn' og imellem enkel-citationstegn, her 'vdpol'. Den femte linie åbner et vindue, hvori de numerisk beregnet løsningskurver plottes (se afsnittet om grafik for videre detaljer). Bemærk, at hver beregningslinie, både her og i M-filen, slutter med et semikolon, hvilket betyder, at de beregnede værdier ikke løbende udskrives på skærmen. Dette er meget tidsbesparende. Ønskes nogle eller alle data og beregninger udskrevet løbende på skærmen undlades de tilsvarende semikolon tegn.

Hvis vi istedet for Van der Pol ligningen ønsker at simulerer følgende modificerede system:

$$\dot{x_1} = x_1(\alpha-x_2^2)-x_2$$

$$\dot{x_2} = \beta x_1$$

$$\dot{x_3} = \gamma \sin(t)$$

for en række forskellige værdier af parameterne $\alpha$, $\beta$ og $\gamma$, skal disse parameter erklæres som værende globale på følgende måde. (Se kapitel 3 om M-filer for mere information om gloable og lokale variable.) M-filen modificeres til:

      function xdot  =  vdpol(t,x)
      global ALPHA BETA GAMMA
      xdot           =  zeros(3,1);
      xdot(1)        =  x(1) .* (ALPHA - x(2).^2)-x(2);
      xdot(2)        =  BETA * x(1);
      xdot(3)        =  GAMMA * sin(t);

Parameterne ALPHA, BETA og GAMMA kan nu ændres interaktivt, f.eks. skrives

 
   tspan    = [0 20];      %Starttid og sluttid 
                           %for integrationen
   xstart   = [0;0.25;0];  %startbetingelse
   global ALPHA BETA GAMMA %parameterne erklæres globale
   ALPHA    = 1.5;         %ALPHA tildeles en værdi
   BETA     = 0.85;        %BETA tildeles en værdi
   GAMMA    = 0;           %GAMMA tildeles en værdi
   [t,x]    = ode23('vdpol',tstart,tslut,xstart);
                           %simuleringen udføres
   plot(t,x)               %kurvene x(1) og x(2) 
                           %plottes versus t

Bemærk, at den modificerede M-fil stadig, noget misvisende, kaldes vdpol.m Herefter kan simuleringer for andre værdier af ALPHA, BETA og GAMMA køres ved blot at tildele disse nye værdier interaktivt efterfulgt af en ny integration og af plot funktionen.

Ønskes et plot af f.eks. x(1) erstattes sidste linie ovenfor blot af

   plot(t,x(:,1));        %kurven x(1) plottes versus t 

Bemærk, at den returnerede variabel x er en matrix med tre søjler og ligeså mange rækker som antallet af skridt i den aktuelle numeriske integration (hvilket ofte er et relativt stort antal, her 185). Ønsker man at kende antallet af rækker skrives blot whos, hvorved oplysninger om bl.a. x matricens dimension udskrives. Når der i ovenstående plot funktion står (t,x(:,1)) som argument, er inputtet til plot funktionen t-vektoren og x(1)-vektoren. Disse har samme dimension og plot funktionen plotter da punkterne (t(1),x(1,1)),(t(2),x(2,1)),..., (t(m),x(m,1)), hvor m er antallet af rækker i x matricen. Da punkterne ligger meget tæt opnås en illusion om at de danner en glat kurve.

Ønskes et faseplot af f.eks. x(2) versus x(1) erstattes sidste linie med

    plot(x(:,1),x(:,2))     %faseplot af x(2) versus x(1)

$\Box$

Det skal nævnes at her blot er medtaget nogle få af de muligheder MATLAB byder på med hensyn til simulering af ordinære (ikke lineære) differentialliningssystemer, for videre brug henvises til manualen.