Modellbildung und Simulation, SS 2010
Blatt 4: ODE-Baukasten

> restart;

> with(DEtools):

> with(LinearAlgebra):

> with(plots):

>

Matrix zu gegebenen Eigenwerten und -vektoren

r: Liste aus den beiden Eigenvektoren

lambda : Liste aus den beiden Eigenwerten

> mkA := proc(r,lambda)
 local T;

 T := <r[1] | r[2]>:

 T . DiagonalMatrix(<lambda[1],lambda[2]>) . MatrixInverse(T)

end;

mkA := proc (r, lambda) local T; T := `<|>`(r[1], r[2]); T . LinearAlgebra:-LinearAlgebra(`<,>`(lambda[1], lambda[2])) . LinearAlgebra:-LinearAlgebra(T) end procmkA := proc (r, lambda) local T; T := `<|>`(r[1], r[2]); T . LinearAlgebra:-LinearAlgebra(`<,>`(lambda[1], lambda[2])) . LinearAlgebra:-LinearAlgebra(T) end procmkA := proc (r, lambda) local T; T := `<|>`(r[1], r[2]); T . LinearAlgebra:-LinearAlgebra(`<,>`(lambda[1], lambda[2])) . LinearAlgebra:-LinearAlgebra(T) end procmkA := proc (r, lambda) local T; T := `<|>`(r[1], r[2]); T . LinearAlgebra:-LinearAlgebra(`<,>`(lambda[1], lambda[2])) . LinearAlgebra:-LinearAlgebra(T) end procmkA := proc (r, lambda) local T; T := `<|>`(r[1], r[2]); T . LinearAlgebra:-LinearAlgebra(`<,>`(lambda[1], lambda[2])) . LinearAlgebra:-LinearAlgebra(T) end proc

> Test1 := mkA([<1,0>,<0,1>],[2,3]);

Test1 := Matrix([[2, 0], [0, 3]])

> Test2 := mkA([<1,1>,<-1,1>],[2,3]);

Test2 := Matrix([[5/2, (-1)/2], [(-1)/2, 5/2]])

> Eigenvectors(Test2);

Vector[column]([[2], [3]]), Matrix([[1, -1], [1, 1]])

>

DGL-System zu gegebenem Richtungsfeld

Neuer Parameter ist der Gleichgewichtspunkt gg.

Solange x nur ein Name ist, sind x[1] und x[2] gute Namen f?r die abh?ngigen Variablen...

> mkDgl := proc(r,lambda,gg)
 local A1, A2;

 A1 := mkA(r,lambda):

 A2 := A1 . (<x[1](t), x[2](t)> - gg):

 {seq(diff(x[i](t),t)=A2[i],i=1..2)}

end;

mkDgl := proc (r, lambda, gg) local A1, A2; A1 := mkA(r, lambda); A2 := A1 . `<,>`(x[1](t), x[2](t))-gg . ; {seq(diff(x[i](t), t) = A2[i], i = 1 .. 2)} end procmkDgl := proc (r, lambda, gg) local A1, A2; A1 := mkA(r, lambda); A2 := A1 . `<,>`(x[1](t), x[2](t))-gg . ; {seq(diff(x[i](t), t) = A2[i], i = 1 .. 2)} end procmkDgl := proc (r, lambda, gg) local A1, A2; A1 := mkA(r, lambda); A2 := A1 . `<,>`(x[1](t), x[2](t))-gg . ; {seq(diff(x[i](t), t) = A2[i], i = 1 .. 2)} end procmkDgl := proc (r, lambda, gg) local A1, A2; A1 := mkA(r, lambda); A2 := A1 . `<,>`(x[1](t), x[2](t))-gg . ; {seq(diff(x[i](t), t) = A2[i], i = 1 .. 2)} end procmkDgl := proc (r, lambda, gg) local A1, A2; A1 := mkA(r, lambda); A2 := A1 . `<,>`(x[1](t), x[2](t))-gg . ; {seq(diff(x[i](t), t) = A2[i], i = 1 .. 2)} end procmkDgl := proc (r, lambda, gg) local A1, A2; A1 := mkA(r, lambda); A2 := A1 . `<,>`(x[1](t), x[2](t))-gg . ; {seq(diff(x[i](t), t) = A2[i], i = 1 .. 2)} end proc

> mkDgl([<1,0>,<0,1>],[2,3],<0,0>);

{diff(x[1](t), t) = 2*x[1](t), diff(x[2](t), t) = 3*x[2](t)}

> mkDgl([<1,1>,<-1,1>],[2,3],<0,0>);

{diff(x[1](t), t) = 5/2*x[1](t)-1/2*x[2](t), diff(x[2](t), t) = -1/2*x[1](t)+5/2*x[2](t)}

>

Beispiele

Eigenvektoren und Eigenwerte:

> r := [<1,0>,<0,1>];

> lambda := [-1,-2];

r := [Vector[column]([[1], [0]]), Vector[column]([[0], [1]])]

lambda := [-1, -2]

r := [<1,0>,<0,1>];
lambda := [-1,2];

r := [<1,0>,<0,1>];
lambda := [1,2];

Gleichgewichtspunkt und Anfangswert:

> gg := <1,1>;

> x0 := gg + 2*r[1] + 4*r[2];

> anfbed := [[x[1](0)=x0[1],x[2](0)=x0[2]]];

gg := Vector[column]([[1], [1]])

x0 := Vector[column]([[3], [5]])

anfbed := [[x[1](0) = 3, x[2](0) = 5]]

Das ist das zugeh?rige DGL-System...

> s := mkDgl(r,lambda,gg);

s := {diff(x[1](t), t) = -x[1](t)+1, diff(x[2](t), t) = -2*x[2](t)+2}

... und seine Koeffizientenmatrix.

> mkA(r,lambda);

Matrix([[-1, 0], [0, -2]])

Plotbereich:

> prange := x[1]=0..5, x[2]=0..5, dirgrid=[11,11];

prange := x[1] = 0 .. 5, x[2] = 0 .. 5, dirgrid = [11, 11]

L?sungskurve und Richtungsfeld:

> lsg := DEplot(s, [x[1],x[2]], t=0..10, anfbed, prange, arrows=NONE,linecolor=BLUE, stepsize=0.01):

> feld := dfieldplot(s,[x[1],x[2]],t=0..1,prange,arrows=MEDIUM, color=BLACK):

> display([lsg,feld]);

[Plot]

Und die Komponenten der L?sung, aufgetragen gegen t:

> p1 := DEplot(s, [x[1],x[2]], t=0..10, anfbed,  linecolor=BLUE, arrows=NONE, stepsize=0.01, scene=[t,x[1]]):

> p2 := DEplot(s, [x[1],x[2]], t=0..10, anfbed, linecolor=RED, arrows=NONE, stepsize=0.01, scene=[t,x[2]]):

> display([p1,p2]);

[Plot]

Variablen l?schen:

> r,lambda,gg,x0 := 'r','lambda','gg','x0';

r, lambda, gg, x0 := r, lambda, gg, x0

>

>