Układy równań różniczkowych

Działania na macierzach

Wprowadzanie macierzy jako zmiennych

Tak jak mówiliśmy na pierwszych zajęciach - macierze to listy list.

m={{(Subscript[a, 1,1]),(Subscript[a, 1,2])},{(Subscript[a, 2,1]),(Subscript[a, 2,2])}}

MatrixForm[m]

m//MatrixForm

Macierze można tworzyć również automatycznie.

m=Table[(Subscript[a, i,j]),{i,1,5},{j,1,3}];
m//MatrixForm

m=Array[(Subscript[a, #1,#2])&,{5,3}];
m//MatrixForm

Aby dostać się do wnętrza listy (macierzy) można użyć notacji [[ ]].

m[[2]]

m[[All,3]]//MatrixForm

m[[2,1]]

m[[2;;4,1;;2]]//MatrixForm

Do sprawdzenia wymiarów macierzy - funkcja Dimensions

Dimensions[m]

Podstawowe działania

A={{-1, 4, -2},{6, 2, -10}};
B={{2, -4, 8},{7,4,2}};

3A-9B//MatrixForm

A={{0, 4, 5},{-5, -1, 5}};
B={{-3, 4}, {-5,-4},{1,-4}};

A.B//MatrixForm

Transpose[A]//MatrixForm

A={{-2,-1,1},{2,1,0},{3,1,-1}};

Det[A]

Minors[A]

Inverse[A]//MatrixForm

Macierze specjalne

Mathematica oferuje całą gamę macierzy “specjalnych” - macierz jednostkową, macierz diagonalną.

IdentityMatrix[5]//MatrixForm

DiagonalMatrix[{2,5,3,1}]//MatrixForm

ConstantArray[2,{3,2}]//MatrixForm

Wartości i wektory własne

Mathematica bez problemu znajduje wartości własne, wektory własne (czy też sam wielomian charakterystyczny).

A={{4,-6},{3,-7}};

eigval=Eigenvalues[A]

eigvec=Eigenvectors[A]

A.eigvec[[2]]-eigval[[2]]*eigvec[[2]]

charpol=CharacteristicPolynomial[A,λ]

Factor[charpol]

Można od razu całość informacji uzyskać jedną funkcją.

Eigensystem[A]

Jeśli chcemy znać jedynie kierunek wektorów i znormalizować długość:

Norm[eigvec[[1]]]

Normalize[eigvec[[1]]]

“Zwykłe” układy równań

Układy równań liniowych można rozwiązywać zwykłą funkcją Solve, ale są funkcje przeznaczone specjalnie do równań liniowych:

A={{-1,1,2},{2,2,0},{0,3,-3}};
b={4,2,0};
x=LinearSolve[A,b]

A.x

A do rozwiązywania równań jednorodnych

A={{3,5,2},{4,7,5},{1,1,-4}};

NullSpace[A]

Macierze jako funkcje

Elementy macierzy nie muszą być liczbami, ale także funkcjami - operacje analityczne są na nich wykonywane element po elemencie.

A={{Cos[3t], Sin[3t], Exp[-t]},{t,t Sin[(t^2)], Exp[t]}};
A//MatrixForm

D[A,t]//MatrixForm

Integrate[A,t]//MatrixForm

Rozwiązywanie układów równań różniczkowych

Równania wyższych rzędów jako układu równań niższych rzędów

Mamy równanie x’’+(x^2-1)x’+x = 0, x(0) = 1, x’(0) = 1, zapiszmy je jako układ równań:
x’ = y
y’ = x’’ = (1-x^2)x’-x = (1-x^2)y-x
z warunkami początkowymi x(0) = 1, y(0) = 1.

Clear[x];

f1[x_,y_]:=y
f2[x_,y_]:=(1-(x^2))y-x

sol=NDSolve[{x'[t]==f1[x[t],y[t]],y'[t]==f2[x[t],y[t]],x[0]==1,y[0]==1},{x[t],y[t]},{t,0,15}]

solutions=Evaluate[{x[t],y[t]}/.sol];
Plot[solutions,{t,0,15},PlotLegends->{"x(t)","y(t)"}]

Show[
StreamPlot[{f1[x,y],f2[x,y]},{x,-5,5},{y,-5,5}],
ParametricPlot[solutions,{t,0,15},PlotStyle->Red]
]

Postać macierzowa układu

Każdy układ równań liniowych można zapisać w postaci macierzowej.
x’’+6x’+5x = cos(t)
Wtedy
x’ = y
y’ = x’’ = -6x’-5x+cos(t)
czyli
X’ = AX + B
gdzie X = (x,y), A = ((0, 1), (-5, -6)), B = (0, cos(t))

Zastosujemy nieco inną składnię w funkcji DSolve - użyjemy tzw. “pure functions”.

A={{0,1},{-5,-6}};
B={0,Cos[t]};
sol=DSolve[{x'[t],y'[t]}==A.{x[t],y[t]}+B,{x,y},t];
FullSimplify[sol]

Dzięki użyciu tych funkcji “czystych” mamy możliwość różniczkowania funkcji, które są rozwiązaniami.

solutionsX=Table[Evaluate[x'[t]/.sol/.{C[1]-> i,C[2]-> 0}],{i,-2,2}];
solutionsY=Table[Evaluate[y[t]/.sol/.{C[1]-> i,C[2]-> 0}],{i,-2,2}];

GraphicsGrid[
{{Plot[solutionsX,{t,0,8}],Plot[solutionsY,{t,0,8}]}}
]

Układ fundamentalny rozwiązań

Jak wyznaczyć układ fundamentalny rozwiązań?

A={{5,-1},{0,3}};

eigval=Eigenvalues[A]
eigvec=Eigenvectors[A]

Z faktu, że macierz A ma dwie rzeczywiste różne wartości własne wynika, że układ fundamentalny ma postać:

Φ=Transpose[{eigvec[[1]] * Exp[eigval[[1]]*t],eigvec[[2]] * Exp[eigval[[2]]*t]}];
Φ//MatrixForm

Φ.{C[1],C[2]}

DSolve[{x'[t],y'[t]}==A.{x[t],y[t]},{x[t],y[t]},t]

Czy to się zgadza? ;-)

Inny przykład, gdy wartości własne są zespolone:

A={{3,-2},{4,-1}};
eigval=Eigenvalues[A]
eigvec=Eigenvectors[A]

Co wtedy? ;-)

Jak rozwiązać zagadnienie początkowe mając macierz fundamentalną?
X(t) = Φ(t) · C
czyli
X(0) = Φ(0) · C ⟹ C = Φ(0)^(-1) · X(0)

A={{1,1},{4,-2}};
eigval=Eigenvalues[A]
eigvec=Eigenvectors[A]

Φ=Transpose[{eigvec[[1]] * Exp[eigval[[1]]*t],eigvec[[2]] * Exp[eigval[[2]]*t]}];
Φ//MatrixForm

c=Inverse[Φ/.t-> 0].{1,-2}

Φ.c

Expand[DSolve[{{x'[t],y'[t]}==A.{x[t],y[t]},x[0]==1,y[0]==-2},{x[t],y[t]},t]]

Funkcja eksponencjalna macierzy

Układ fundamentalny można znaleźć (dla układów równań liniowych) również inaczej.
Wiemy, że dla układów X’ = AX układ fundamentalny to exp(At). Wiemy jednak też, że “zwykła” funkcja Exp działa element po elemencie, a “prawdziwa” wartość exp(At) to suma szeregu
exp(At) = I + At + 1/(2!) A^2 t^2 + 1/(3!) A^3 t^3 + ...

Clear[A]

Series[Exp[A t],{t,0,4}]

W środowisku Mathematica mamy odpowiednią funkcję na to

A={{2,5},{-4,-2}};

Φ=MatrixExp[A t]

czyli rozwiązanie ogólne ma postać

Φ.{C[1],C[2]}

DSolve[{x'[t],y'[t]}==A.{x[t],y[t]},{x[t],y[t]},t]

Jak podnosić macierze do potęgi?

MatrixPower[A,5]