7. Vaja

Diferencialne enačbe prvega reda

Diferencialna enačba je enačba, v kateri nastopajo neodvisna spremenljivka x, funkcija neodvisne spremenljivke y in njeni odvodi. Če v enačbi nastopa samo prvi odvod, je to diferencialna enačba 1. reda, če pa v enačbi nastopajo tudi višji odvodi, pa govorimo o diferencialnih enačbah drugega reda, tretjega reda... Medtem ko pri običajnih enačbah iščemo neznano spremenljivko, pri diferencialnih enačbah iščemo neznano funkcijo. Oglejmo si najprej geometrijski pomen diferencialnih enačb.

0. Naj bo

f[x_,y_]:=1+x^2+(2 x y/(1+x^2));

Radi bi rešili diferencialno enačbo y'=f(x,y) z metodo smernic.
Najprej narišemo nivojski diagram funkcije f(x,y):

cp=ContourPlot[f[x,y],{x,-3,3},{y,-3,3},Contours->Table[i,{i,-10,10}],ContourShading->False]

Nato nivojnice opremimo s kratkimi črticami z ustreznim kotom.

vp=VectorPlot[{1,f[x,y]},{x,-3,3},{y,-3,3},VectorScale->{Automatic,Automatic,None}];
vcp=Show[cp,vp]

Na koncu smernice povežemo v krivulje. Te krivulje so rešitve diferencialne enačbe y'=f(x,y)

DSolve[y'[x]==f[x,y[x]],y[x],x];
Table[(y[x]/.%)/.C[1]->i,{i,-2,2}]//Flatten;
Plot[%,{x,-3,3},PlotRange->{-3,3},PlotStyle->Thickness[0.015]];
Show[vcp,%]

Naš rezultat primerjajmo z rezultatom, ki ga dobimo z že vgrajenim ukazom StreamPlot.

StreamPlot[{1,f[x,y]},{x,-3,3},{y,-3,3},StreamPoints->Coarse]

Najpreprostejše diferencialne enačbe prvega reda so diferencialne enačbe z ločljivimi spremenljivkami. 
To so enačbe oblike y' = f(x,y), kjer lahko funkcijo f(x,y) zapišemo v obliki f(x,y)=g(x) h(y) za neki funkciji g in h.

1. Reši diferencialno enačbo y' sin(x)= y ln(y).

Najprej damo izraz, v katerem nastopa samo y, na levo stran enačbe, izraz, v katerem nastopa samo x, pa na desno stran enačbe. Nato integriramo vsako stran posebej.

leva=Integrate[1/(y Log[y]),y]

desna=Integrate[1/Sin[x],x]+D

Da izrazimo y kot funkcijo x, uporabimo ukaz Solve

Solve[leva == desna, y]

V Mathematici rešujemo diferencialne enačbe z ukazom DSolve.

??DSolve

DSolve[y'[x]Sin[x]==y[x]Log[y[x]],y[x],x]

2. Reši diferencialne enačbe z ločljivimi spremenljivkami s pomočjo ukaza DSolve:
a) y' = sin(x),
b) y' = -(y/x),
c) y' = (x^2)y,
d) y' = (y^2)sin(x).

DSolve[y'[x]== Sin[x],y[x],x]

DSolve[y'[x]== -y[x]/x,y[x],x]

DSolve[y'[x]==x^2 y[x],y[x],x]

DSolve[y'[x]==y[x]^2 Sin[x],y[x],x]

Linearna diferencialna enačba 1. reda je enačba oblike y' = p(x) y + q(x). Če je q(x) = 0, pravimo, da je enačba homogena, drugače pa je nehomogena.

3. Reši linearno diferencialno enačbo (1+(x^2)) y' -2 x y =(1+(x^2))2.

Najprej poiščemo splošno rešitev homogene enačbe, ki je oblike (Subscript[y, h])= c  (e^(∫p(x)dx)).

Clear[x,y,c,p,q]

p[x_]:=(2x)/(1+(x^2))

q[x_]:=1+(x^2)

yh=c  Exp[Integrate[p[x],x]]

Nato poiščemo partikularno rešitev nehomogene enačbe, ki je oblike (Subscript[y, p])= c(x) (e^(∫p(x)dx)), kjer je c(x) =  ∫q(x) (e^(-∫p(x)dx))dx. 
To formulo lahko poenostavimo, če upoštevamo definicijo (Subscript[y, h]). Dobimo (Subscript[y, p])=(Subscript[y, h])∫((q(x))/(Subscript[y, h]))   dx.

yp=Integrate[q[x] / yh,x] yh

Splošna rešitev nehomogene linearne diferencialne enačbe je oblike y = (Subscript[y, h])+(Subscript[y, p]). Dobimo:

y=yh+yp

Poiščimo še rešitev dane diferencialne enačbe s pomočjo ukaza DSolve.

Clear[x,y,c,p,q]

DSolve[(1+x^2)y'[x] - 2 x y[x] ==(1+x^2)^2,y[x],x]

Splošna rešitev diferencialne enačbe prvega reda ima en prosti parameter. V Mathematici so ti parametri označeni s konstantami C[1]. Vsaka izbira konstante C[1] nam določa eno rešitev diferencialne enačbe. Vrednost konstante C[1] lahko določimo s pomočjo začetnih pogojev.

4. V isti koordinatni sistem narišite rešitve diferencialne enačbe y' = -(y/(x+1)), ki zadoščajo začetnim pogojem y(1) = 1, y(1) = 2 in y(1) = 3.

Clear[x,y]

DSolve[y'[x]==-y[x]/(x+1),y[x],x]

Splošna rešitev dane diferencialne enačbe je torej y(x)=((C[1])/(1+x)). Vrednost konstante C[1] pri začetnem pogoju y(1) = 1 določimo iz enačbe 1=y(1)=((C[1])/(1+1)).

Solve[((C[1])/(1+1))==1,C[1]]

Rešitev dane diferencialne enačbe pri začetnem pogoju y(1)=1 je torej y(x)=(2/(1+x)).

Z ukazom  DSolve lahko rešimo diferencialno enačbo tudi pri danih začetnih pogojih.

r1=DSolve[{y'[x]==-y[x]/(1+x),y[1]==1},y[x],x]

r2=DSolve[{y'[x]==-y[x]/(1+x),y[1]==2},y[x],x]

r3=DSolve[{y'[x]==-y[x]/(1+x),y[1]==3},y[x],x]

Plot[{Evaluate[y[x]/.r1],Evaluate[y[x]/.r2],Evaluate[y[x]/.r3]},{x,-5,5}, PlotRange->{-6,6}]

Logistična diferencialna enačba N'=kN(Nmax-N) se uporablja za modeliranje spreminjanja populacije. Pri tem je N trenutno število elementov v populaciji, k koeficient rasti populacije in Nmax maksimalno število elementov v populaciji.

5. Poišči rešitev logistične diferencialne enačbe v odvisnosti od začetne populacije, koeficienta rasti in maksimalne populacije.

DSolve[{y'[x]== k y[x](Nmax-y[x]),y[0]==y0},y[x],x]

Poglejmo si graf rasti populacije pri vrednosti Nmax = 100.

Nmax=100;

Manipulate[Plot[(((E^(k Nmax x)) Nmax y0)/(Nmax-y0+(E^(k Nmax x)) y0)),{x,0,2},PlotRange->{0,100}],{{k,0.1,"Koeficient rasti"},-1,1},{{y0,10,"Začetna populacija"},0,100}]

Za mnoge diferencialne enačbe ne znamo poiskati splošne rešitve. Takrat si pomagamo z metodami za približno določitev tiste rešitve diferencialne enačbe, ki gre skozi določeno točko. Ogledali si bomo Eulerjevo metodo.

Iščemo rešitev y=y(x) diferencialne enačbe y' = f(x,y), ki gre skozi točko ((Subscript[x, 0]),(Subscript[y, 0])). Rešitev iščemo na intervalu [(Subscript[x, 0]),(Subscript[x, 0])+d]. Najprej si izberemo število korakov n in dolžino koraka h=(d/n). Sedaj definirajmo (Subscript[x, 1]) = (Subscript[x, 0])+h in (Subscript[y, 1]) = (Subscript[y, 0])+f((Subscript[x, 0]),(Subscript[y, 0]))h. Točka ((Subscript[x, 1]),(Subscript[y, 1])) leži na premici skozi točko ((Subscript[x, 0]),(Subscript[y, 0])) s predpisanim naklonom y' = f((Subscript[x, 0]),(Subscript[y, 0])), njena x-koordinata pa se od x-koordinate prvotne točke razlikuje za h. Postopek induktivno nadaljujemo s predpisom (Subscript[x, i+1]) = (Subscript[x, i])+h,  (Subscript[y, i+1]) = (Subscript[y, i])+f((Subscript[x, i]),(Subscript[y, i]))h. za i<n. Ko te točke zaporedoma povežemo z daljicami, dobimo lomljeno črto, ki je približek za rešitev diferencialne enačbe. Približek je tem boljši, čim večje je število korakov.

6. Dana je diferencialna enačba y' = (((x^2)+(y^2))/10). Z Eulerjevo metodo konstruiraj približno rešitev enačbe, ki gre skozi točko (1,1) na intervalu [1,5]. Primerjaj rezultate, ki jih dobiš, če izvedeš 3, 5, 10 ali 100 korakov Eulerjeve metode.

Clear[x,y,crte]
f[x_,y_]:=(x^2+y^2)/10
x0=1;
y0=1;
d=4;

Približek dobimo tako, da najprej s pomočjo For zanke definiramo seznam točk ((Subscript[x, i]),(Subscript[y, i])), nato pa te točke povežemo z lomljeno črto.

n=3;
h=d/n;

crte={{x0,y0}}//N;
For[i=1,i<=n,i++,{x,y}=Last[crte];AppendTo[crte,{x+h,y+f[x,y]h}]];
priblizek3=ListPlot[crte,PlotJoined->True]

n=5;
h=d/n;

crte={{x0,y0}}//N;
For[i=1,i<=n,i++,{x,y}=Last[crte];AppendTo[crte,{x+h,y+f[x,y]h}]];
priblizek5=ListPlot[crte,PlotJoined->True]

n=10;
h=d/n;

crte={{x0,y0}}//N;
For[i=1,i<=n,i++,{x,y}=Last[crte];AppendTo[crte,{x+h,y+f[x,y]h}]];
priblizek10=ListPlot[crte,PlotJoined->True]

n=100;
h=d/n;

crte={{x0,y0}}//N;
For[i=1,i<=n,i++,{x,y}=Last[crte];AppendTo[crte,{x+h,y+f[x,y]h}]];
priblizek100=ListPlot[crte,PlotJoined->True]

Poiščimo še natančno rešitev dane diferencialne enačbe z ukazom DSolve in jo narišimo ter primerjajmo približno in točno rešitev.

Clear[res,x,y];
res=DSolve[{y'[x]==(x^2+y[x]^2)/10,y[1]==1},y[x],x];

resitev=Plot[Evaluate[y[x]/.res],{x,1,5},PlotStyle->{Red}]

Show[{priblizek3,priblizek5,priblizek10,priblizek100,resitev}]

Nekaterih diferencialnih enačb Mathematica ne zna rešiti simbolično. V tem primeru lahko numerični približek rešitve diferencialne enačbe, ki ustreza danim začetnim pogojem, poiščemo z ukazom NDSolve.

7. Numerično reši diferencialno enačbo y' = (x^2)-(y^2) pri začetnem pogoju y(0)=0 na intervalu [0,3].

r1=NDSolve[{y'[x]==x^2-y[x]^2,y[0]==0},y,{x,0,3}]

Plot[Evaluate[y[x]/.r1],{x,0,3},PlotRange->All]

Diferencialne enačbe drugega reda

Splošno rešitev diferncialne enačbe drugega reda znamo najti le v zelo posebnih primerih. Podrobno si bomo ogledali linearne diferencialne enačbe drugega reda.

Linearna diferencialna enačba drugega reda je enačba oblike y'' = p(x) y' + q(x) y+r(x). Če sta p in q konstantni funkciji in je r(x) = 0, govorimo o homogeni linearni enačbi drugega reda s kontantnimi koeficienti.

8. Poišči splošno rešitev homogene linearne diferencialne enačbe drugega reda s konstantnimi koeficienti y'' - 4y '+ 3y = 0.

Tako enačbo rešujemo z nastavkom y=(e^(λx)),kjer je λ konstanta,ki jo moramo še določiti. Izračunamo prvi in drugi odvod nastavka ter to vstavimo v prvotno enačbo. Dobimo karakteristično enačbo oblike
 0=(λ^2)-pλ-q in poiščemo njene rešitve.

Solve[l^2-4l+3==0,l]

Če ima enačba dve različni rešitvi (Subscript[λ, 1 ])in (Subscript[λ, 2, ])sta (e^((Subscript[λ, 1])x))in (e^((Subscript[λ, 2])x)) linearno neodvisni rešitvi dane diferencialne enačbe. Tudi vsaka njuna linearna kombinacija reši dano diferencialno enačbo. Splošna rešitev diferencialne enačbe v našem primeru bo torej oblike (C1e^x)+(C2e^(3x)).

Če ima karakteristična enačba eno samo rešitev λ, sta linearno neodvisni rešitvi oblike (e^(λx)) in x (e^(λx)).

Če ima karakteristična enačba kompleksni rešitvi a + bi in a- bi, je splošna rešitev oblike y = C1 (e^ax)cos(bx) + C2 (e^ax)sin(bx).

9. Poišči splošne rešitve naslednjih homogenih linearnih diferencialnih enačb drugega reda s konstantnimi koeficienti:
a) y''-5y'+6y = 0,
b) y''-2y'+y = 0,
c) y''+y = 0,
d) y''-2y'+10y = 0.

Poglejmo si zgoraj omenjene možnosti neodvisnih rešitev na danih primerih s pomočjo ukaza  DSolve.

DSolve[y''[x]-5y'[x]+6y[x]==0,y[x],x]

DSolve[y''[x]-2y'[x]+y[x]==0,y[x],x]

DSolve[y''[x]+y[x]==0,y[x],x]

DSolve[y''[x]-2y'[x]+10y[x]==0,y[x],x]

10. Poišči splošno rešitev diferencialne enačbe y'' -2 y' + 5 y = (e^x)

Splošna rešitev nehomogene diferencialne enačbe je vsota splošne rešitve homogene enačbe in partikularne rešitve nehomogene enačbe. Partikularno rešitev dobimo iz splošne rešitve homogene enačbe z variacijo konstante. Če je y = C1 y1 + C2 y2 splošna rešitev homogene enačbe, velja:
C1' y1+C2' y2=0;
C1' y1'+C2' y2'=r;
Odtod izračunamo, da je 
C1' = -((r (Subscript[y, 2]))/((Subscript[y, 1]) (Subscript[y, 2])'-(Subscript[y, 1])' (Subscript[y, 2]))) in C2' = ((r (Subscript[y, 1]))/((Subscript[y, 1]) (Subscript[y, 2])'-(Subscript[y, 1])' (Subscript[y, 2]))).   
C1(x) in C2(x) dobimo z integriranjem C1' in C2'.
Partikularna rešitev je (Subscript[y, p]) =C1(x)(Subscript[y, 1])(x) + C2(x) (Subscript[y, 2]) (x).

Clear[yh,yp]

Solve[l^2-2 l +5 ==0,l]

f[x_,C1_,C2_]:=C1 Exp[x] Cos[2x]+C2 Exp[x] Sin[2x]

yh=f[x,C1,C2]

y1[x_]:=  Exp[x] Cos[2x]
y2[x_]:=  Exp[x] Sin[2x]

D1=Integrate[Exp[x] y2[x]/(y1[x]  y2'[x]-y1'[x] y2[x]),x]
D2=Integrate[Exp[x] y1[x]/(y1[x]  y2'[x]-y1'[x] y2[x]),x]

yp=f[x,D1,D2]

y=yp+yh

Rešimo dano diferencialno enačbo še s pomočjo ukaza  DSolve.

Clear[y]

DSolve[y''[x]-2y'[x]+5y[x]==E^x,y[x],x]

11. Poišči splošno rešitev nehomogenih diferencialnih enačb:
a) y'' -5 y' + 6 y = 1,
b) y'' -5 y' + 6 y = (e^x),
a) y'' -5 y' + 6 y = sin(x).

DSolve[y''[x]-5y'[x]+6y[x]==1,y[x],x]

DSolve[y''[x]-5y'[x]+6y[x]==E^x,y[x],x]

DSolve[y''[x]-5y'[x]+6y[x]==Sin[x],y[x],x]

12. Enačba dušenega nihanja mx''(t) + cx'(t) +kx(t) = 0 je osnoven primer homogene diferencialne enačbe drugega reda s konstantnimi koeficienti. Pri tem je m masa uteži, c koeficient dušenja in k koeficient vzmeti. Poglejmo si primer, ko vzmet raztegnemo do začetne amplitude A in jo natu spustimo.

DSolve[{ m y''[x]+c y'[x]+k y[x]==0,y[0]==A,y'[0]==0},y,x]

f[x_,A_,c_,k_,m_]:=((A (-c (E^(((-c-(Sqrt[(c^2)-4 k m])) x)/(2 m)))+c (E^(((-c+(Sqrt[(c^2)-4 k m])) x)/(2 m)))+(E^(((-c-(Sqrt[(c^2)-4 k m])) x)/(2 m))) (Sqrt[(c^2)-4 k m])+(E^(((-c+(Sqrt[(c^2)-4 k m])) x)/(2 m))) (Sqrt[(c^2)-4 k m])))/(2 (Sqrt[(c^2)-4 k m])))

Animacijo gibanja vzmeti si lahko ponazorimo z naslednjim programom K. F. Stephensa, ki je dostopen na portalu MathWorld.

Manipulate[
Grid[{{
(* sinusoidal plot *)
Plot[-UnitStep[t-τ](1+f[τ,A,c,k,m]),{τ,0,30},PlotStyle->{ColorData["HTML","SlateBlue"],Thickness[.005]},PlotRange->{0,-2},ImageSize->{450,350},Frame->True,Axes->False,FrameTicks->{Automatic,{{0,+1},{-1,0},{-2,-1}}}],
(* mass and spring graphics *)
Graphics[{ParametricPlot[{0.5+0.1Sin[8x],-(x/(2π))(1+f[t,A,c,k,m])},{x,0,2π},PlotStyle->Directive[Thick,GrayLevel[.2]],Axes->None,AspectRatio->1,ImageSize->{112,350}][[1]],ColorData["HTML","SlateBlue"],Rectangle[{0.3,-0.9-f[t,A,c,k,m]},{.7,-1.1-f[t,A,c,k,m]}]},Frame->True,PlotRange->{0,-2},ImageSize->{112,280},FrameTicks->{False,{{0,+1},{-1,0},{-2,-1}}}]}}],
{{A,.3,"Začetni odmik"},0.1,1,Appearance->"Labeled"},
{{k,2,"Togost vzmeti"},1,5,Appearance->"Labeled"},
{{m,2,"Masa"},1,5,Appearance->"Labeled"},{{c,0.2,"Koeficient dušenja"},0,1,Appearance->"Labeled"},{{t,(10^(-6)),"Čas"},(10^(-6)),30,ControlType->Trigger,DefaultDuration->(10^(-6)),DisplayAllSteps->True,AnimationRate->2},AutorunSequencing->{{4,30}}]

13. Kot zgled nehomogene diferencialne enačbe drugega reda si poglejmo poševni met v homogenem polju sile teže. Točka med padanjem opiše parabolično krivuljo, ki jo podamo v parametrični obliki (x(t),y(t)). Obliko funkcij x in y dobimo iz drugega Newtonovega zakona, ki v tem primeru pravi x''(t)=0 in y''(t)=-g, kjer je g težni pospešek. Začetni pogoji so x(0)=x0, x'(0)=vx, y(0)=y0,y'(0)=vy.

DSolve[{x''[t]==0,x[0]==x0,x'[0]==vx},x[t],t]

DSolve[{y''[t]==-g,y[0]==y0,y'[0]==vy},y[t],t]

Poglejmo si še program Roba Morrisa in Oleksandra Pavlyka, ki grafično demonstrira poševni met.

Clear[h]

Manipulate[
Module[{x,y,soln,eqns,t,xmax,ymax},

eqns={x''[t]==0,x'[0]==v*Cos[θ],x[0]==0,y''[t]==-g,y'[0]==v*Sin[θ],y[0]==h};
soln=Flatten[Quiet@NDSolve[eqns,{x,y},{t,0,Infinity},Method->{"EventLocator","Event"->y[t],"EventAction":>Throw[tf=t,"StopIntegration"],"Direction"->-1},MaxSteps->Infinity]];

ymax=First@Quiet@FindMaximum[{y[t]/.soln,0≤t≤tf},{t,tf/2}];
xmax=x[tf]/.soln;

Labeled[Show[{ParametricPlot[{x[t],y[t]}/.soln,{t,0,Min[tf,end]},AxesOrigin->{0,0},ImageSize->{400,300},PlotRange->{{0,Max[1.1xmax,1.1ymax]},{0,Max[1.1xmax,1.1ymax]}},AxesLabel->{Style["x(t)",Italic,Bold,"Label"],Style["y(t)",Italic,Bold,"Label"]},ImagePadding->30],Graphics[{Darker[Green],Disk[(({x[Min[tf,end]],y[Min[tf,end]]}/.soln)),Scaled[.01]]}]}],Style[Row[{"Čas padanja: ",tf," s",Spacer[30],"Maksimalna višina: ",ymax," m",Spacer[30],"Domet: ",xmax," m"}],"Label"]]],{{v,10,"Začetna hitrost"},.1,20,Appearance->"Labeled"},{{θ,Pi/4.,"Začetni kot"},.001,Pi/2.,Appearance->"Labeled"},{{h,10,"Začetna višina"},.001,20,Appearance->"Labeled"},{{g,9.8,"Gravitacijski pospešek"},3,20,Appearance->"Labeled"},{{end,.001,"Prikaži"},.001,Dynamic[tf],ControlType->Trigger},{{end,.001},.001,2.49,ControlType->None},AutorunSequencing->{7},SaveDefinitions->True,Initialization:>{tf=Infinity},TrackedSymbols:>{θ,h,g,end}]

14. Reši z Eulerjevo metodo diferencialno enačbo y'' = y na intervalu [0,5] pri začetnih pogojih y(0)=1 in y'(0)=1. Rešitev primerjaj z eksaktno rešitvijo in pa z rešitvijo ki jo dobiš s pomočjo ukaza NDSolve.

Tudi diferencialne enačbe drugega reda lahko rešujemo z Eulerjevo metodo.

Iščemo rešitev y=y(x) diferencialne enačbe y'' = f(x,y,y'), pri začetnih pogojih y((Subscript[x, 0]))=(Subscript[y, 0]) in y'((Subscript[x, 0]))=(Subscript[z, 0]). Rešitev spet iščemo v n korakih na intervalu [(Subscript[x, 0]),(Subscript[x, 0])+d] ; dolžina enega koraka je h=(d/n). Induktivno definirajmo (Subscript[x, i+1]) = (Subscript[x, i]) +h in (Subscript[y, i+2]) = 2(Subscript[y, i+1]) -(Subscript[y, i]) +(h^2)f ((Subscript[x, i]), (Subscript[y, i]),  (((Subscript[y, i+1])-(Subscript[y, i]))/h)) za i<n-1. Ko te točke zaporedoma povežemo z daljicami, dobimo lomljeno črto, ki je približek za rešitev dane diferencialne enačbe.

Clear[x,y,crte,f];
f[x_,y_,z_]:=y;
x0=0;
y0=1;
d=5;

Izračunali bomo približek za rešitev, ki ga dobimo s pomočjo 10 korakov.

n=10;
h=d/n;
x1=x0+h;
y1=y0+h;

crte={{x0,y0},{x1,y1}}//N;
For[i=0,i<=n-2,i++,{{x,y}=crte[[i+1]],z=Last[crte[[i+2]]]};AppendTo[crte,{x+2h,2z-y+f[x,y,(z-y)/h]h^2}]];
priblizek10=ListPlot[crte,PlotJoined->True,AxesOrigin->{0,0}]

Eksaktna rešitev dane diferencialne enačbe je funkcija y(x)=(e^x).

eksaktnar=Plot[E^x,{x,0,5},PlotStyle->{Red}]

Izračunajmo še numerično rešitev enačbe s pomočjo ukaza NDSolve in vse tri rešitve narišimo na en graf.

Clear[x,y]

r=NDSolve[{y''[x]==y[x],y[0]==1,y'[0]==1},y[x],{x,0,5}]

numr=Plot[Evaluate[y[x]/.r],{x,0,5},PlotStyle->{Green},PlotRange->All]

Show[{priblizek10,eksaktnar,numr}]

Sistem diferencialnih enačb prvega reda

V uporabi pogosto naletimo tudi na sisteme diferencialnih enačb. Sistem dveh diferencialih enačb prvega reda ima obliko y' = f(x,y,z), z' = g(x,y,z).

Linearna diferencialna enačba drugega reda je enačba oblike y'' = p(x) y' + q(x) y+r(x). Če sta p in q konstantni funkciji in je r(x) = 0, govorimo o homogeni linearni enačbi drugega reda s kontantnimi koeficienti.

11. Poišči splošno rešitev sistema diferencialnih enačb y' = -z, z' = y.

Sistem rešimo tako, da iz prve enačbe izrazimo z = -y' in vstavimo v drugo enačbo. Dobimo enačbo -y'' = y. Karakteristiki te enačbe sta +i in -i, zato je splošna rešitev y(x) = C1cos(x)+C2sin(x).

DSolve[-y''[x]==y[x],y[x],x]

Sedaj lahko izračunamo z = -y'. Torej je z(x) = C1sin(x)-C2cos(x).

Clear[x,z]

y[x_]:=C[1] Cos[x]+C[2] Sin[x]

z[x_]=-y'[x]

Rešimo sistem še s pomočjo ukaza DSolve.

Clear[x,y,z]

DSolve[{y'[x]==-z[x],z'[x]==y[x]},{y[x],z[x]},x]

Naloge

1. Reši na dva načina diferencialni enačbi:
a) y' = (y^3),
b) y' = y+x.

2. Reši diferencialno enačbo y' = (x^3)+(y^2) z Eulerjevo metodo na intervalu [0,1] pri začetnem pogoju y(0)=1. Primerjaj rezultate, ki jih dobiš, če izvedeš 3, 5 in 10 korakov Eulerjeve metode z rešitvijo, ki jo dobiš s pomočjo ukaza DSolve.

3. Numerično reši diferencialno enačbo y' = sin((y^2))x pri začetnem pogoju y(0)=1 na intervalu [-3,3].

4. a) Poišči splošno rešitev diferencialne enačbe 4y'' + 4y' + y = 4+x,    
b) Poišči rešitev diferencialne enačbe y'' + 4y' + 4y = 0, ki zadošča začetnemu pogoju y(0)=2, y'(0)=0.

5. Reši sistem diferencialnih enačb prvega reda y' = y + z, z' = y - z.