MatLab - løsning av system av differensialligninger

[krje1980]

Hei.

Jeg holder på å lære meg MatLab på egen hånd. Har et spørsmål som jeg håper noen MatLab-flinkiser kan svare meg på! I tekstboken jeg bruker er det en egen seksjon om hvordan man kan løse ordinære differensialligninger. Jeg har forstått det meste her, men ser ikke helt logikken i et eksempel boken bruker for å illustrere bruken av såkalt Event Location function. Dvs at man skriver et program hvor man instruerer om at dersom en viss begivenhet inntreffer så skal f.eks. MatLab avslutte utregningen av differensialligningene selv om vi opprinnelig har spesifisert et større domene.

Eksempelet i boken begynner med at man definerer en funksjon som følger:

Kode: Velg alt

function [value, isterminal, dircn] = proj(t,z,flag);
    g = 9.81;
    if nargin < 3 | isempty(flag)
        value = [z(2); 0; z(4); -g];
    else
        switch flag
            case 'events'
                value = z(3);
                isterminal = 1;
                dircn = -1;
            otherwise
                error('function not programmed for this event');
        end
    end

Deretter går man frem som følger:

Kode: Velg alt

tspan = [0 6];
z0 = [0, 5*cos(pi/4), 0, sin(pi/4)];
options = odeset('events','on');
[t, y] = ode45('proj',tspan,z0,options)

MatLab vil nå skrive fem kolonner med tall, hvorav den første kolonnen representerer verdier for t, mens de andre kolonnene representerer verdiene for z(1), z(2), z(3) og z(4) ved korresponderende t. Når z(3) treffer verdien 0, avsluttes utregningene.

Dersom jeg ikke hadde tatt med "options" linjen ovenfor, og kun skrevet:

Kode: Velg alt

[t, y] = ode45('proj',tspan,z0)

Så får jeg ut akkurat de samme kolonnene som over, men nå fortsetter utregningene etter at z(3) = 0, og helt til vi når endepunktet for tspan.

Det jeg imiderltd ikke skjønner er hvordan MatLab kan regne ut verdiene for z(1), z(2), z(3) og z(4) når vi aktiverer "events". Slik jeg tolker logikken i funksjonen "proj", så vil jo man da, når man aktiverer options under utførelsen av ode45, forholde seg til det som står i "else"-delen av funksjonen. Men det er jo under "if"-delen at selve systemet beskrives under

Kode: Velg alt

value = [z(2); 0; z(4); -g];

Hvorfor er det da slik at MatLab løser systemet helt fint selv om man aktiverer "events"?

Dersom noen kan forklare meg logikken her så vil jeg være meget takknemlig! Dette har plaget meg i hele natt!

[drgz]

Jeg får at isempty(flag) slår til med begge kodene du har lagt fram på første iterasjon - og da vil du jo alltid gå inn i if-delen av koden til å begynne med. Men etter første iterasjon så går den over til else-delen av koden noen ganger - så da må det være enten pga flag = 'events'.

Iterasjon 1: flag = '', uansett om du har med odeopts eller ikke
Iterasjon 2: flag = 'events'
osv

Du kan sette en debug-marker f.eks på linje 2 i proj.m og kjøre koden. Så er det bare å steppe med F10 helt til koden er ferdig for å se hva som skjer. Hvis du f.eks på første iterasjon setter flag = 'events' manuelt mens du debugger, så vil du se at koden ikke klarer å fullføre ganske fort - men utover det tror jeg du må se litt etter selv da jeg ikke har brukt ode45 noe særlig før.

[krje1980]

drgz:

Tusen takk for svar! Kan det være slik at dette er noe som rett og slett ligger innebygd i ode45 metoden? At den først etter å ha gjennomført en iterasjon, aktiverer "events" for å sjekke om "eventen" har funnet sted? Og så gjentas dette for hver verdi for t? Det vil i så fall kunne forklare hvorfor programmet oppfører seg som det gjør - at "eventen" aktiveres først etter at programmet har kjørt gjennom if-delen.

[drgz]

Er usikker på hvorfor. Prøvde å se om events = 'on' som default, men det ser ikke ut til det. Ellers ser jeg at hvis du f.eks setter options = odeset('events','off') så kjører koden til t = 6, men det har, utfra hva jeg kan se, ingenting å gjøre med hva som skjer på første iterasjonen gjennom proj.m - der er flag = '' uansett hva du gjør.

Og nå så jeg at du skrev at koden går til t = 6 når du tar vekk linjen med options også . Tydeligvis på tide å legge seg her:p

[krje1980]

OK. Setter stor takk på hjelpen. Jeg er ganske sikker på at dette må ha noe med selve aktiveringen av ode45 å gjøre. Det er jo åpenbart at if-delen gjennomføres for hver verdi i tspan frem til "eventen" inntreffer. Samtidig merker jeg at når først "eventen" inntreffer (at z(3) = 0), så vil også denne linjen med vektorverdier skrives ut i vinduet før iterasjonene avsluttes. Igjen kan dette tyde på at det først er ETTER verdiene er kalkulert at programmet sjekker om "eventen" har funnet sted.