Differentiaaliyhtälöiden ratkaiseminen ja havainnollistaminen

1 Vektorikenttä

Harjoitustehtävissä 1.3 ja 1.4 tarkasteltavaa vektorikenttää N voi havainnollistaa VectorPlot - komennolla

VectorPlot[{-x/Sqrt[x^2+y^2],-y/Sqrt[x^2+y^2]},{x,-1,1},{y,-1,1}]

Automaattisesti valitut pisteet, joissa vektorit piirretään, eivät välttämättä ole tarkasteltavan vektorikentän kannalta paras valinta . Tällöin voidaan käyttää VectorPoints-optiota. Seuraavassa valitaan pisteet r ympyrältä kulman 2π/M välein.

M=24;r=4;
Pisteet = Flatten[Table[N[{m Cos[k 2 Pi/M]/r, m  Sin[k 2 Pi/M]/r}], {k, 0, M}, {m, 1, r}], 1];
VectorPlot[{-x/Sqrt[x^2+y^2],-y/Sqrt[x^2+y^2]},{x,-1,1},{y,-1,1},VectorPoints->Pisteet]

2 Tason differentiaaliyhtälö

Tarkastellaan tason differentiaaliyhtälöä (OverDot[x])=Bx. Etsitään ensin ominaisarvot ja ominaisvektorit. Eigensystem
antaa ominaisarvot 2 ja 1 ja vastaavat ominaisvektorit {3,1} ja {2,1}.

B={{3,5},{1,-1}}; es=Eigensystem[B]

KK=Transpose[es[[2]]]

Inverse[KK].B.KK

Piirretään vektorikenttä neliössä [-1, 1]*[-1, 1]

v=VectorPlot[{3x+5y,x-y},{x,-1,1},{y,-1,1},VectorColorFunction->None,VectorStyle->Red]

Piirretään ominaisavaruuksien osa, joka on neliössä [-1,1]×[-1,1]

e=Graphics[{Line[{{1,1/5},{-1,-1/5}}],Line[{{-1,1},{1,-1}}]}]

Piirretään 25 radan osia  neliössä [-1, 1]*[-1, 1]

r=StreamPlot[{3x+5y,x-y},{x,-1,1},{y,-1,1},StreamPoints->25]

Yhdistetään edellä piirretyt havainnollistukset yhteen kuvaan

Show[{e,r,v}]

Esimerkki spiraalinielu

K={{3,1},{1,1}};
A=K.{{-1,2},{-2,-1}}.Inverse[K]

v=VectorPlot[{-5x+10y,-2x+3y},{x,-1,1},{y,-1,1},VectorColorFunction->None,VectorStyle->Red]

v2=VectorPlot[{-5x+10y,-2x+3y},{x,-2,2},{y,-1,1},VectorPoints->Map[K.#&,Pisteet],AspectRatio->Automatic,VectorMarkers->"Dart",VectorSizes->.5]

s=StreamPlot[{-5x+10y,-2x+3y},{x,-2,2},{y,-2,2},StreamPoints->30]

Show[{v2,s}]

Eigensystem[A]

M={{2,-1},{1,0}}

Inverse[M].A.M

Esimerkki 3.4: Vaimeneva värähtelijä

Dynaaminen moduli antaa mahdollisuuden tarkastella vaimenevan värähtelijän ratkaisuja parametrien k ja b eri arvoilla . Kappaleen massa m on kiinnitetty mutta sitäkin voi vaihtaa koodissa. Parametritasoon on piirretty osa paraabelia, jolla systeemi on kriittisesti vaimennettu.

m=1;back=Plot[2 Sqrt[x],{x,0,5},AspectRatio->1];

DynamicModule[{k={1,0}},
Column[{
Row[{LocatorPane[Dynamic[k],back],Slider2D[Dynamic[k],{0,5},Appearance->"Labeled",ImageSize->Medium]}],
Row[{Dynamic[StreamPlot[{y,-k[[1]]x+-k[[2]]y},{x,-2,2},{y,-2,2},StreamPoints->30]],"     ",  
Column[{"{k,b}="Dynamic[k],
Row[{"Ominaisarvot ",Dynamic[Eigenvalues[{{0,1},{-k[[1]]/m,-k[[2]]/m}}]]
}]
}]
}]
}]
]

Jos emme välitä fysikaalisesta uskottavuudesta voimme sallia myös negatiivisen jousuvakion ja kitkan .

m=1;back=ParametricPlot[{Abs[x],Sign[x]2 Sqrt[Abs[x]]},{x,-5,5},PlotRange->{{-5,5},{-5,5}},AspectRatio->1];

DynamicModule[{k={1,0}},
Column[{
Row[{LocatorPane[Dynamic[k],back],Slider2D[Dynamic[k],{0-5,5},Appearance->"Labeled",ImageSize->Medium]}],
Row[{Dynamic[StreamPlot[{y,-k[[1]]x+-k[[2]]y},{x,-2,2},{y,-2,2},StreamPoints->30]],"     ",  
Column[{"{k,b}="Dynamic[k],
Row[{"Ominaisarvot ",Dynamic[Eigenvalues[{{0,1},{-k[[1]]/m,-k[[2]]/m}}]]
}]
}]
}]
}]
]

Harjoitustehtävä 3.7

DynamicModule[{k=0},
Column[{
Row[{Slider[Dynamic[k],{0-5,5},Appearance->"Labeled",ImageSize->Medium]}],
Row[{Dynamic[StreamPlot[{k x+2y,-2x},{x,-2,2},{y,-2,2},StreamPoints->30]],"     ",  
Column[{"α="Dynamic[k],
Row[{"Ominaisarvot ",Dynamic[Eigenvalues[{{k,2},{-2,0}}]]
}]
}]
}]
}]
]

Harjoitus 3.8

b=-1/4;k=1;m=1;
dy7[x0_,y0_]:=NDSolve[{x'[t]==y[t],y'[t]==-k x[t]/m- b y[t]/m,x[0]==x0,y[0]==y0},{x[t],y[t]},{t,-10,10}];
ratk=dy7[1,0];
pss=ParametricPlot[Evaluate[{x[t],y[t]}/.ratk],{t,0,10}, PlotRange->{{-4,3},{-3,2}}, PlotPoints->100,AxesLabel->{"(Subscript[x, 1])","(Subscript[x, 2])"},PlotStyle->{Black,Thickness[0.005]}];

b=-4;k=1;m=1;
dy7[x0_,y0_]:=NDSolve[{x'[t]==y[t],y'[t]==-k x[t]/m- b y[t]/m,x[0]==x0,y[0]==y0},{x[t],y[t]},{t,0,20}];
ratk=dy7[1,0];
pc=ParametricPlot[Evaluate[{x[t],y[t]}/.ratk],{t,0,1}, PlotRange->{{-6,3},{-30,2}}, PlotPoints->100,AxesLabel->{"(Subscript[x, 1])","(Subscript[x, 2])"},PlotStyle->{Red,Thickness[0.005]}];

b=-6;k=1;m=1;
dy7[x0_,y0_]:=NDSolve[{x'[t]==y[t],y'[t]==-k x[t]/m- b y[t]/m,x[0]==x0,y[0]==y0},{x[t],y[t]},{t,0,1}];
ratk=dy7[1,0];
ps=ParametricPlot[Evaluate[{x[t],y[t]}/.ratk],{t,0,1}, PlotRange->{{-4,3},{-30,2}}, PlotPoints->100,AxesLabel->{"(Subscript[x, 1])","(Subscript[x, 2])"},PlotStyle->{Blue,Thickness[0.005]}];

Show[{pss,pc,ps},PlotRange->{{-4,3},{-10,2}}]

5 Epälineaariset differentiaaliyhtälöt

StreamPlot[{-x/2+y+(x^2+y^2)x/2,-x-y/2+(x^2+y^2)y/2},{x,-1.5,1.5},{y,-1.5,1.5},StreamPoints->{{{{1,0},Red},{{.9,0},Green},{{.8,0},Red}, 500}},StreamColorFunction->None]