Toets PiMM vrijdag 18 maart 2016, 15-17.15u

Notebooks voor je begint eerst copiëren naar je desktop. Onderweg notebooks regelmatig saven ivm mogelijke crash. Als je klaar bent eerst saven naar je desktop, daarna pas de gesavede notebooks submitten via Blackboard. Check zelf dat de submit gelukt is.

Vragen: Vragen kun je alleen per email stellen aan de docent I.H.M.vanStokkum@uva.nl Die beslist of hij eventueel een hint geeft, waarmee dan later in de beoordeling rekening wordt gehouden.

Communicatie : Het is niet toegestaan om via het internet te communiceren, of te zoeken op het internet.

Hints: Lees de opgaven eerst goed door, maak een plan, en begin daarna pas met het intypen van commando's. Controleer je werk kritisch.

Werk nauwkeurig. Let op veel voorkomende fouten zoals:

teveel of te weinig spaties en accenten (alleen een spatie gebruiken voor vermenigvuldiging), 
verwisseling van letter o en cijfer 0, of verwisseling van letter l en cijfer 1
functies beginnen met een hoofdletter (bv Sin, NDSolve, PlotStyle etc)
argumenten van functies staan tussen rechte haken (blokhaken), en worden gescheiden door komma’s
door te clicken tussen twee cellen en dan te typen krijg je automatisch een nieuwe input cell (NB probeer nooit input te typen in een output cell )
vergeten van tijdsargument, of verkeerd tijdsargument bij een beginvoorwaarde,
vergeten termen van een noemer bij elkaar te nemen dmv ronde haken,
verwisseling van plus en min tekens, 
het verschil tussen = (Set) en == (Equal) en := (SetDelayed, voor definitie van een eigen functie met argumenten) 
het verschil tussen /. en //. (ReplaceRepeated)
slordig gebruik van globale variabelen 
vergeten van ; (CompoundExpression) in Module
zorg dat je { } goed matchen, en ook je [ ] en je ( )

Mocht je onverhoopt in een niet-eindigende berekening terecht komen, selecteer dan uit het Evaluation menu Abort Evaluation ( Alt+. ).

Gebruik indien nodig lijstindicering [[ ]] (Part) of functies als Flatten, First, Last om resultaten "zonder accolades" te krijgen. Alleen in uiterste noodgevallen kun je ook Copy en Paste gebruiken.

in noodgevallen: selecteer dan uit het Evaluation menu Quit Kernel
of begin met een schone lei  (verwijderen van alle globale variabelen):

Remove["Global`*"]

Als dit statement bovenaan je notebook staat kun je om alle globale variabelen te verwijderen ook Evaluate Notebook selecteren uit het Evaluation menu

Wanneer Manipulate crasht ...

Evaluation / Quit Kernel
Evaluation / Dynamic Updating Enabled

Hints

Wanneer je onverhoopt ergens vastloopt, merk dan op dat sommige deelopgaven ook zonder sommige tussenresultaten gemaakt kunnen worden. Wanneer je problemen hebt met kleuren, plot dan zonder. Wanneer je de syntax van een commando wel weet, maar er is onderweg iets mis gegaan, geef dan toch die syntax.

Waardering staat bij elk onderdeel na de letter aangegeven.Totaal te behalen punten 100.  Afronding gebeurt pas op het eind van het gehele vak Practicum Natuurkunde 1 B.

Model voor een biologische oscillator

Remove["Global`*"]

(a4) Definieer een functie g met vier argumenten genaamd u, v, j en k. g is een breuk met als teller 
2uk
en als noemer 
 v-u+vJ+uK+(Sqrt[(((v-u+vj+uk))^2)-4(v-u)uk])

g[u_,v_,j_,k_]:=2u k/(v-u+v j+u k+Sqrt[(v-u+v j+u k)^2-4(v-u) u k])

(b2) maak een object ratesR gelijk aan 
k0 maal g(R,1,k1,k1) plus S minus R minus R maal X 
(NB hier geen tijdsafhankelijkheid)

ratesR=k0 g[R,1,k1,k1]+S-R-R X

(c2) maak een object ratesX gelijk aan 
k2 maal R minus k3 maal X

ratesX=k2 R-k3 X

(d1) maak een nieuw object ratesRt waarbij je in ratesR een tijdsafhankelijkheid toevoegt voor R en X (dwz R vervangt door R(t) en X vervangt door X(t))

ratesRt=ratesR/.{R->R[t],X->X[t]}

(e1) maak een nieuw object ratesXt waarbij je in ratesX een tijdsafhankelijkheid toevoegt voor R en X

ratesXt=ratesX/.{R->R[t],X->X[t]}

(f1) Maak een lijst valuesf met substitutieregels met waarden k0=4, k1=0.3, k2=0.1, k3=0.075, S=0.2

valuesf={k0->4,k1->0.3,k2->0.1,k3->0.075,S->0.2}

(g12) Los de gekoppelde differentiaalvergelijkingen R' = ratesRt en X' = ratesXt numeriek op met beginvoorwaarden R (0) = 0.8 en X (0) = 1 voor t tot 100.

rule=NDSolve[{R'[t]==ratesRt,X'[t]==ratesXt,R[0]==0.8,X[0]==1}/.valuesf,{R,X},{t,0,100}]//Flatten

(h7) Maak een plot van R en X als functie van de tijd. Kleur de curven blauw en rood, geef ze een Thickness van 0.01, en gebruik de optie AxesLabel.

{rt,xt}={R[t],X[t]}/.rule

Plot[{rt,xt},{t,0,100},PlotStyle->{{Blue,Thickness[0.01]},{Red,Thickness[0.01]}},AxesLabel->{"t","rt,xt"}]

(2 punten voor objecten, 2 voor plot, 3 voor opties)

(i5) Maak een parametrische plot van R en X als functie van de tijd. Kleur de curve zwart, geef  een Thickness van 0.01, en gebruik de optie AxesLabel (met R en X langs de juiste as). Gebruik voor al dit soort phase plots (faseplaatjes) als PlotRange voor R van 0 tot 3, en voor X van 0 tot 1.5.

ploti=ParametricPlot[{rt,xt},{t,0,100},PlotStyle->{Black,Thickness[0.01]},AxesLabel->{"R","X"},PlotRange->{{0,3},{0,1.5}}]

(2 puntenvoor plot, 3 voor opties)

(j10) Schrijf een functie twoplot met argumenten R0 en X0 voor de beginvoorwaarden en tend voor de eindtijd van de berekening (en van de plaatjes). Gebruik in die functie Module om de beide voorgaande plaatjes naast elkaar te tekenen.

twoplot[R0_,X0_,tend_]:=Module[{}, g[u_,v_,j_,k_]:=2u k/(v-u+v j+u k+Sqrt[(v-u+v j+u k)^2-4(v-u) u k]);
 ratesR=k0 g[R,1,k1,k1]+S-R-R X;ratesX=k2 R-k3 X;ratesRt=ratesR/.{R->R[t],X->X[t]};ratesXt=ratesX/.{R->R[t],X->X[t]};
 values={k0->4,k1->0.3,k2->0.1,k3->0.075,S->0.2}; rule=NDSolve[{R'[t]==ratesRt,X'[t]==ratesXt,R[0]==R0,X[0]==X0}/.values,{R,X},{t,0,tend}]//Flatten;
 {rt,xt}={R[t],X[t]}/.rule;ploth=Plot[{rt,xt},{t,0,tend},PlotStyle->{{Blue,Thickness[0.01]},{Red,Thickness[0.01]}},AxesLabel->{"t","rt,xt"}];
ploti=ParametricPlot[{rt,xt},{t,0,tend},PlotStyle->{Black,Thickness[0.01]},AxesLabel->{"R","X"},PlotRange->{{0,3},{0,1.5}}];{ploth,ploti}]

Test je functie  twoplot met argumenten 1, 0.8 en 100, en met argumenten 1.8, 0.2 en 150.

twoplot[0.8,1,100]

twoplot[1.8,0.2,150]

(k3) Breid je functie twoplot uit tot een nieuwe functie twoplotS met een argument S0 voor de stimulus sterkte S. zorg mbv de juiste substitutie dat dit argument S0 gebruikt wordt om de oplossing te berekenen.

twoplotS[R0_,X0_,tend_,S0_]:=Module[{}, g[u_,v_,j_,k_]:=2u k/(v-u+v j+u k+Sqrt[(v-u+v j+u k)^2-4(v-u) u k]);
 ratesR=k0 g[R,1,k1,k1]+S-R-R X;ratesX=k2 R-k3 X;ratesRt=ratesR/.{R->R[t],X->X[t]};ratesXt=ratesX/.{R->R[t],X->X[t]};
 values={k0->4,k1->0.3,k2->0.1,k3->0.075,S->S0}; rule=NDSolve[{R'[t]==ratesRt,X'[t]==ratesXt,R[0]==R0,X[0]==X0}/.values,{R,X},{t,0,tend}]//Flatten;
 {rt,xt}={R[t],X[t]}/.rule;ploth=Plot[{rt,xt},{t,0,tend},PlotStyle->{{Blue,Thickness[0.01]},{Red,Thickness[0.01]}},AxesLabel->{"t","rt,xt"}];
ploti=ParametricPlot[{rt,xt},{t,0,tend},PlotStyle->{Black,Thickness[0.01]},AxesLabel->{"R","X"},PlotRange->{{0,3},{0,1.5}}];{ploth,ploti}]

(l 12) Onderzoek nu dit model beschreven in twoplotS mbv Manipulate als functie van X0 (van 0 tot 1.5), R0 (van 0 tot 3), S0 (van 0 tot 0.5) en tend (van 100 tot 300).

Manipulate[twoplotS[R0,X0,tend,S0],{R0,0,3},{X0,0,1.5},{S0,0,0.5},{tend,100,300}]

Tussen welke grenswaarden (op een honderdste nauwkeurig) van S0 vind je oscillaties ? kies een waarde van S0 waar je oscillaties vindt, en varieer nu de beginvoorwaarden R0 en X0. wat constateer je? beschrijf nu in je eigen woorden wat een limietcyclus is.

tussen 0.066 en 0.407 (0.07 en 0.41 is ook goed) bij alle beginvoorwaarden R0 en X0 vind je de oscillaties. een limietcyclus is dus een periodieke beweging waar dit systeem naar toe wordt getrokken wanneer S0 tussen 0.07 en 0.41 ligt), onafhankelijk vd beginvoorwaarden,

Isoclines

(m6) bereken vanwege de gecompliceerde R afhankelijkheid van ratesR de beide isoclines als functie van X.

RisoX=Solve[ratesR==0,X]//Flatten

XisoX=Solve[ratesX==0,X]//Flatten

(n6) Maak een plot (genaamd plotn) van de beide isoclines (waar ratesR en ratesX nul zijn) waarbij je R varieert van 0 tot 3, gebruik valuesf en als PlotRange voor de beide isoclines van 0 tot 1.5. Kleur de curven  rood en blauw, geef ze een Thickness van 0.01, en gebruik de optie AxesLabel.

valuesf={k0->4,k1->0.3,k2->0.1,k3->0.075,S->0.2}

plotn1= Plot[Evaluate[{X/.RisoX,X/.XisoX}/.valuesf],{R,0,3},PlotRange->{0,1.5},PlotStyle->{{Blue,Thickness[0.01]},{Red,Thickness[0.01]}},AxesLabel->{"R","X"}]

laat plotn en de parametrische plot gemaakt mbv twoplot[2.7, 0.1,50] samen zien.

plotn2=twoplot[2.7,0.1,50]

Show[plotn1,plotn2[[2]]]

(o6) Schrijf een functie isoplot met argument S0 waarin je de beide isoclines tekent als functie van de stimulus sterkte S. Gebruik in die functie Module. test isoplot met S0 gelijk aan 0.01 en 0.49.

isoplot[S0_]:=Module[{},ratesR=k0 g[R,1,k1,k1]+S-R-R X;ratesX=k2 R-k3 X;RisoX=Solve[ratesR==0,X]//Flatten; XisoX=Solve[ratesX==0,X]//Flatten;values={k0->4,k1->0.3,k2->0.1,k3->0.075,S->S0};plotn1= Plot[Evaluate[{X/.RisoX,X/.XisoX}/.values],{R,0,3},PlotRange->{0,1.5},PlotStyle->{{Blue,Thickness[0.01]},{Red,Thickness[0.01]}},AxesLabel->{"R","X"}]]

isoplot[0.01]

isoplot[0.49]

(p5) Schrijf nu een nieuwe functie twoplotSiso waarbij je de isoplot toevoegt aan de parametrische plot van twoplotS. test met de argumenten  0.737412,0.983216,100,0.2

twoplotSiso[R0_,X0_,tend_,S0_]:=Module[{}, g[u_,v_,j_,k_]:=2u k/(v-u+v j+u k+Sqrt[(v-u+v j+u k)^2-4(v-u) u k]);
 ratesR=k0 g[R,1,k1,k1]+S-R-R X;ratesX=k2 R-k3 X;ratesRt=ratesR/.{R->R[t],X->X[t]};ratesXt=ratesX/.{R->R[t],X->X[t]};
 values={k0->4,k1->0.3,k2->0.1,k3->0.075,S->S0}; rule=NDSolve[{R'[t]==ratesRt,X'[t]==ratesXt,R[0]==R0,X[0]==X0}/.values,{R,X},{t,0,tend}]//Flatten;
 {rt,xt}={R[t],X[t]}/.rule;ploth=Plot[{rt,xt},{t,0,tend},PlotStyle->{{Blue,Thickness[0.01]},{Red,Thickness[0.01]}},AxesLabel->{"t","rt,xt"}];
ploti=ParametricPlot[{rt,xt},{t,0,tend},PlotStyle->{Black,Thickness[0.01]},AxesLabel->{"R","X"},PlotRange->{{0,3},{0,1.5}}];RisoX=Solve[ratesR==0,X]//Flatten; XisoX=Solve[ratesX==0,X]//Flatten;plotn1= Plot[Evaluate[{X/.RisoX,X/.XisoX}/.values],{R,0,3},PlotRange->{0,1.5},PlotStyle->{{Blue,Thickness[0.01]},{Red,Thickness[0.01]}},AxesLabel->{"R","X"}];{ploth,Show[ploti,plotn1]}]

twoplotSiso[0.737412,0.983216,100,0.2]

(q6) Bereken het stationaire punt van deze gekoppelde differentiaalvergelijkingen (waarbij je valuesf gebruikt) en geef de coordinaten van het stationaire punt (met positieve R en X waarden) als startwaarden voor de passende aanroep van twoplotSiso met tend 120.

statrules=Solve[{ratesR==0,ratesX==0}/.valuesf,{R,X}]

{Rs,Xs}={R,X}/.statrules[[2]]

twoplotSiso[Rs,Xs,120,0.2]

(r5) wat neem je waar? wat concludeer je over de aard van het stationaire punt?

ongeveer 70 s blijft het systeem vlak in de buurt van het stationaire punt, en daarna vliegt het naar de limietcyclus. het stationaire punt is instabiel. dit komt door de eindige nauwkeurigheid van zestien decimalen.

save nu eerst je notebook. 
om de correctheid en stand alone bruikbaarheid van twoplotSiso te testen dien je eerst alle globale variabelen te verwijderen. Doe dit door de volgende text cel mbv Format Style in een input te veranderen (shortcut Alt-9) en uit te voeren.

Remove["Global`*"]

(s6) Herhaal nu je definitie van twoplotSiso hieronder en test door  aanroep van twoplotSiso met argumenten 0.737412,0.983216,100,0.2. Als alles werkt klopt dit, en krijg je weer plaatjes. anders kun je altijd nog je zojuist gesavede notebook inleveren.

twoplotSiso[R0_,X0_,tend_,S0_]:=Module[{}, g[u_,v_,j_,k_]:=2u k/(v-u+v j+u k+Sqrt[(v-u+v j+u k)^2-4(v-u) u k]);
 ratesR=k0 g[R,1,k1,k1]+S-R-R X;ratesX=k2 R-k3 X;ratesRt=ratesR/.{R->R[t],X->X[t]};ratesXt=ratesX/.{R->R[t],X->X[t]};
 values={k0->4,k1->0.3,k2->0.1,k3->0.075,S->S0}; rule=NDSolve[{R'[t]==ratesRt,X'[t]==ratesXt,R[0]==R0,X[0]==X0}/.values,{R,X},{t,0,tend}]//Flatten;
 {rt,xt}={R[t],X[t]}/.rule;ploth=Plot[{rt,xt},{t,0,tend},PlotStyle->{{Blue,Thickness[0.01]},{Red,Thickness[0.01]}},AxesLabel->{"t","rt,xt"}];
ploti=ParametricPlot[{rt,xt},{t,0,tend},PlotStyle->{Black,Thickness[0.01]},AxesLabel->{"R","X"},PlotRange->{{0,3},{0,1.5}}];RisoX=Solve[ratesR==0,X]//Flatten; XisoX=Solve[ratesX==0,X]//Flatten;plotn1= Plot[Evaluate[{X/.RisoX,X/.XisoX}/.values],{R,0,3},PlotRange->{0,1.5},PlotStyle->{{Blue,Thickness[0.01]},{Red,Thickness[0.01]}},AxesLabel->{"R","X"}];{ploth,Show[ploti,plotn1]}]

twoplotSiso[0.737412,0.983216,120,0.2]

beschrijf wat je ziet in dit plaatje, en verklaar het verschil met het plaatje in onderdeel q

ongeveer 20 s blijft het systeem vlak in de buurt van het stationaire punt, en daarna vliegt het naar de limietcyclus (op een andere manier). dit komt omdat de beginvoorwaarden maar 6 decimalen nauwkeurig gelijk zijn aan die van het stationaire punt, terwijl ze in onderdeel q zestien decimalen nauwkeurig waren.