Analytisk løsning af differentialligninger

MuPad

De seneste versioner af MATLAB kommer med et Computer Algebra System (CAS) kaldet MuPAD (Tidligere versioner af MATLAB kommer i stedet med Maple), som startes ved at skrive,
>> pad = mupad
i MATLABs kommandovindue.

Hvis man vil sende information mellem MuPAD og MATLAB er det nødvendigt at man kalder mupad med en variabel, i dette tilfælde ``pad''. MuPAD kan bruges til symbolsk beregning og grafisk opsætning af artikler og afleveringer, den store fordel er dog integrationen med MATLAB.

Vi vil i de første opgaver undersøge raketligningen analytisk, start derfor MuPAD, og læs ``Getting Started'' sektionen i den indbyggede hjælp. MuPAD har to input formater, et til tekst og billeder, og et til beregninger. Beregningsinput har en lille firkantet kasse foran sig. Når man tager udtryk ind i beregningsområdet kan det evalueres ved at trykke [enter].

Lad os definere en funktion \(f\) og plotte den i et interval \(x=1...10\)

f:=k*exp(2/x)
plot(f|k=4,x=1..10)
Bemærk at ingen argumenter er specificeret for \(f\), det gør det muligt løbende at bestemme hvad der er variable og hvad der er konstanter. I plot funktionen specificeres konstanterne med en en substitutionssteg: | efter funktionen og variablens interval specificeres med \(x=start..slut\). Se figur figuren ovenfor. \newpage Vi kan løse ligning~\ref{rocketeq} i MuPAD med ode::solve() der som navnet antyder løser ordinære differentialligninger. Det kan anbefales at læse vejledningen til ode:: funktionerne da de kan meget mere end hvad der er vist her.
[caption={Løsning af differentialligning (ligning~\ref{rocketeq}) med MuPAD}, label=lst:ode]
f := ode::solve({m(t)*v'(t)=-v_x*m'(t), v(0)=v0},v(t));	
>> {v0 - v_x*ln(m(t)) + v_x*ln(m(0))}
I listing~\ref{lst:ode} er ligning~\ref{rocketeq} specificeret, hvor \(v'(t)\) svarer til \(dv/dt\), \(m(t)\) fortæller at m er afhængig af \(t\) og \(v_x\) er en konstant. Der er også defineret en begyndelsesbetingelse, nemlig at \(v(0)=v_0\). Det sidste argument \(v(t)\) er den funktion vi ønsker at finde. Resultatet svarer nu til ligning~\ref{eq:dv} hvis vi flytter \(v0\) over på venstresiden. Før vi kan arbejde videre med funktionen bør vi først erstatte \(m(t)\) og \(m(0)\), da vi ikke har en tidsafhængigt massefunktion, men bare en start og slut værdi.
[caption={Substitution af udtryk i MuPAD}, label=lst:mupadsub]
f :=f|m(t)=m1|m(0)=m0

I listing~\ref{lst:mupadsub} benytter vi igen | til at beskrive hvad de enkelte led skal erstattes med, og skriver resultatet tilbage til variablen \(f\). Funktionen kan nu retuneres til MATLAB som et symbolsk objekt ved at benytte getVar() kommandoen.

[caption={Import af udtryk fra MuPAD til MATLAB. bemærk at udtrykket skal eksistere på forhånd ved at man har eksekveret MuPAD scriptet før man kalder getVar(). variabelnavnet pad blev defineret i listning~\ref{lst:nypad}.},label=lst:getvar]
>> fm = getVar(pad, 'f')
  fm = v0 + v_x*log(m0) - v_x*log(m1)
>> dv = matlabFunction(f)
  dv = @(m0,m1,v0,v_x)reshape([v0+v_x.*log(m0)-v_x.*log(m1)],[1,1]) 

Den første linje importerer MuPADs variabel f til en symbolsk variabel fm i MATLAB. matlabFunction(f) konverterer et symbolsk MATLAB udtryk til en anonym MATLAB funktion, som f.eks. kan bruges til numeriske beregninger eller plotting (reshape() kommandoen gør i dette tilfælde ikke noget, men er indsat af MATLAB under konveteringen fra symbolsk udtryk).

Opgaver

  1. Løs ligning \ref{rocketeq} analytisk med MuPAD som beskrevet ovenfor.
  2. Importer funktionen fra MuPAD til MATLAB med \lstinline!matlabFunction(.)!.
  3. Benyt ezcontourf() funktionen til at undersøge funktionen for forskelige forhold \(\mu = m_0/m\) og mundingshastigheder \(v_{ex}=0...3000\) m/s. (Hint: du skal definere m som et forhold af m0 for at kunne plotte funktionen i 2 dimensioner: ezcontourf(@(x,y)dv(x,1, 1/y), [1, 2000], [0, 30]))
  4. Løsning ligning~\ref{rocketeq} for \(m(t)\) i stedet for \(v(t)\), med begyndelsesbetingelsen \(m(0)=m0\). (løsning: \(m_1(\Delta v) = m_0 e^{-\Delta v/v_{ex}}\)). Benyt substitutioner og funktionen simplify() i MuPAD til at forme resultatet som ønsket.
  5. Den typiske mundingshastighed på en Falcon 9 raket er \(v_{ex} = 2837\) m/s. Plot funktionen \(m_1(\Delta v)\) for \(v_{ex} = 2837\) m/s enten i MuPAD eller MATLAB. Start massen for en Falcon 9 raket er \(m_0 = 415700\) kg, hvor stor en hastighedsforskel kan man opnå hvis \(90\%\) af raketten er brændstof?