Programozásról

Kevesebb szorzás kell:

Horner-alak

HornerForm[11 (x^3)-4 (x^2)+7 x+2]

Racionálias függvényre is:

HornerForm[(11 (x^3)-4 (x^2)+7 x+2)/((x^2)-3x+1)]

Többváltzósra is

HornerForm[(11 (x^3)y-4 (x^2)+7 x+2(y^2))(y (x^2)-3x+1)]

Paraméteres együtthatókkal is

HornerForm[(11 (x^3)y-4 p (x^2)+7 q x+2(y^2))(y (x^2)-3(p^2)x+r),{x,y}]

Karatsuba-algoritmus

ClearAll[a,b,c,d];
p1=(a x+b)(c x+d)//Expand

4 szorzás

Expand[u=(a+b)(c+d)]

p2=a c (x^2)+(u-a c-b d)x+b d

p1==p2

Simplify@%

Karatsuba, 3 szorzás

Mátrixszorzás

Még ügyesebbek is lehetünk.

Hyperlink["Strassen","http://en.wikipedia.org/wiki/Strassen_algorithm"]

Map - megírjuk

mymap1[f_,l_List]:=Table[f[x],{x,l}]

mymap2[f_,l_]:=mymap2[f,{},l];(*lánckezdés, inicializálás, iniciáció*)
mymap2[g,{a,s,d}]

mymap2[f_,{x___},{y_,z___}]:=mymap2[f,{f[y],x},{z}];(*elongáció, folytatás*)
mymap2[g,{cc,dd},{a,s,d}]

mymap2[_,x_,{}]:=x;(*elongáció, folytatás*)
mymap2[g,{cc,dd,a,s,d}]

Pókháló diagram

Első megoldás: Riffle

A lényeges rész ez:

{{(Subscript[x, 0]),(Subscript[x, 0])},{(Subscript[x, 0]),F[(Subscript[x, 0])]},{F[(Subscript[x, 0])],F[(Subscript[x, 0])]},{F[(Subscript[x, 0])],F[F[(Subscript[x, 0])]]},{F[F[(Subscript[x, 0])]],F[F[(Subscript[x, 0])]]},{F[F[(Subscript[x, 0])]],F[F[F[(Subscript[x, 0])]]]},{F[F[F[(Subscript[x, 0])]]],F[F[F[(Subscript[x, 0])]]]}}

NestList megadja a pályát:

palya=NestList[F,(Subscript[x, 0]),3]

Amit kapni akarunk, az két sorozat összefésüléséből ered, az egyik a szögfelezőn találtható pontokból,

{{(Subscript[x, 0]),(Subscript[x, 0])},{F[(Subscript[x, 0])],F[(Subscript[x, 0])]},{F[F[(Subscript[x, 0])]],F[F[(Subscript[x, 0])]]},{F[F[F[(Subscript[x, 0])]]],F[F[F[(Subscript[x, 0])]]]}}

a mási az F függvény gráfján található pontokból áll:

{{(Subscript[x, 0]),F[(Subscript[x, 0])]},{F[(Subscript[x, 0])],F[F[(Subscript[x, 0])]]},{F[F[(Subscript[x, 0])]],F[F[F[(Subscript[x, 0])]]]}}

Itt van az összeféfülő függvény:

Riffle[{aa,xx,cc},{dd,ee,ff}]

Állítsuk elő az összefésülendőket:

tpp=Transpose[{palya,palya}]

tmr=Transpose[{Most@palya,Rest@palya}]

Riffle[tpp,tmr]

Második megoldás: Partition

Vannak alkalmas opcionális paraméterei: az offszet.

Partition[{qq,ww,ee,rr,tt,zz,uu,ii},2,1]

Partition[palya,2,1]

Avagy

Partition[Flatten[Transpose[{palya,palya,palya}]],2]

Most térünk vissza

7. A súgó alábbi példáját tanulmányozva

nds=NDSolve[{x''[t] + x[t]==0,x'[0]==π,x[0]==0, WhenEvent[{x[t]==2,x[t]==-3},x'[t]->-x'[t]]},x,{t,0,20}];
Plot[{x[t]/.nds,2,-3},{t,0,20}]

modellezzük egy szoba fűtését: amikor a hőmérséklet eléri a 22 fokot, a fűtés kikapcsol, és a levegő a Newton-féle lehűlési törvény szerint hűl. Amint eléri a 18 fokot, bekapcsolunk egy 40 fokos kályhát. Rajzoljunk le 4  periódust.

nd=First@NDSolve[{D[T'[t]==-(T[t]-18),t],T[0]==21,T'[0]==18-T[0],WhenEvent[{T[t]==18.01},T'[t]->40-T[t]],WhenEvent[{T[t]==21.99},T'[t]->18-T[t]]},T,{t,0,30}]

Plot[T[t]/.nd,{t,0,30}]

Oláh Judit

8. Szoba fűtése

nds=NDSolve[{x''[t] + x[t]==0,x'[0]==π,x[0]==0, WhenEvent[{x[t]==2,x[t]==-3},x'[t]->-x'[t]]},x,{t,0,20}];
Plot[{x[t]/.nds,2,-3},{t,0,20}]

Newton lehűlési törvénye
dT/dt = -k*(T-T0)

ww[t_]=DSolveValue[{x'[t]==-40*(x[t]-18), x[0]==22},x[t],t]

D[ww[t],{t,2}]

ww1[t_]=D[ww[t],t]

ww1[0]

ff[t_]=DSolveValue[{x'[t]==0.04*(-x[t]+40), x[0]==18},x[t],t]

D[ff[t],{t,2}]

ff1[t_]=D[ff[t],t]

nds1=NDSolve[{x''[t] ==1600x[t]-28800,x'[0]==-160,x[0]==22, WhenEvent[{x[t]==18, x[t]==22},{x'[t]->-0.04*x[t]+1.6, x'[t]->-40x[t]+720}]},x,{t,0,20}];
Plot[{x[t]/.nds1,18,22},{t,0,2}, Method->{"EquationSimplification"->"Residual"},PlotLabel->"A szoba hőmérséklete", AxesLabel->{"Idő", "Hőmérséklet C"}]

nds2=NDSolve[{x'[t] ==-40x[t]+720,x[0]==22, WhenEvent[{x[t]==18, x[t]==22},{x'[t]->-0.4x[t]+16, x'[t]->-40x[t]+720}]},x,{t,0,20}];
Plot[{x[t]/.nds2,18,22},{t,0,5}, Method->{"EquationSimplification"->"Residual"},PlotLabel->"A szoba hőmérséklete", AxesLabel->{"Idő", "Hőmérséklet C"}]

Hibajelzést ad.
Azt írja, nem állítható x’[t] a Whenevent paranccsal, mert elsőrendűnél nincs magasabb derivált.

nds3=ParametricNDSolve[{x'[t] ==-b*(x[t]-a),x[0]==18, WhenEvent[{x[t]==22,x[t]==18},{{a->40, b->-400},{a->10, b->400}}]},x,{t,0,10}, {a, b},MaxSteps->∞];
Plot[{x[t]/.nds3,18,22},{t,0,10},PlotLabel->"A szoba hőmérséklete", AxesLabel->{"Idő", "Hőmérséklet C"}]

8. Az α és β mátrix ismeretében írjuk föl egy reakció kinetikai differenciálegyenletét tömeghatás típusú kinetikával. Megadom a koordinátás alakot, de szeretnék koordinátamentes megoldást is a haladó vegyészektől.

(Subscript[(OverDot[c]), m])=(∑_(r=1)^R)(β(m,r)-α(m,r))(Subscript[k, r])(∏_(p=1)^M)(Subscript[c, p]^(α(p,r)))

MatrixForm[α={{1,1,0},{0,1,1}}]

MatrixForm[β={{2,0,0},{0,2,0}}]

M=2;R=3;
(Subscript[f, m_]):=Sum[(β[[m,r]]-α[[m,r]])(Subscript[k, r])Product[((Subscript[c, p])^(α[[p,r]])),{p,1,M}],{r,1,R}]

(Subscript[f, 1])

(Subscript[f, 2])

(({(Subscript[c, 1]),(Subscript[c, 2])})^(Transpose[α]))

pt=Power[{(Subscript[c, 1]),(Subscript[c, 2])},#]&/@Transpose[α]

tpt=Times@@@pt

ktpt={k1,k2,k3} tpt

(β-α).ktpt

rhs[a_,b_,k_]:=(b-a).(k Times@@@(Power[{x,y},#]&/@Transpose[a]))

rhs[α,β,{(Subscript[k, 1]),(Subscript[k, 2]),(Subscript[k, 3])}]

de=Thread[{x´[t],y´[t]}==(rhs[α,β,{(Subscript[k, 1]),(Subscript[k, 2]),(Subscript[k, 3])}]/.{x->x[t],y->y[t]})]

Mathematica  és differenciálegyenletek

Tóth János, Analízis Tanszék, jtoth@math.bme.hu, www.math.bme.hu/~jtoth, H 311

Innen http://www.math.bme.hu/~jtoth/Mma/M_M_2008.pdf letölthető a könyvünk: Szili, L., Tóth, J.: Matematika és Mathematica, ELTE Eötvös Kiadó, Budapest,1996.

Itt lehet tanulni a programról: 
BMETE929012	    	Felsőbb “Mathematica”
BMETE925015 	Számítógép az alkalmazott analízisben
BMETE92MC12	Informatika MC

Első félév, második félév, első félév, a harmadik angolul, lassabb, 4 óra, az első kettő magyarul, 2 óra.

Mathematica

Steven Wolfram, részecskefizika, sejtautomaták, Mathematica, A New Kind of Science, wolframalpha.com

Szimbolikus megoldások

Most: köz. de. igen, PDE nem

Nagyjából azt tudja, amit Kamke.

E. Kamke : Differentialgleichungen : Lösungsmethoden und Lösungen, 1. Gewöhnliche Differentialgleichungen, 2. Partielle  Differentialgleichungen

Egyszerű példák

Lineáris egyenletek és rendszerek

Állandó együtthatós inhomogén lineáris egyenlet

ClearAll[egy,kezd1,kezd2];
egy=y'[x]==2x-y[x];
kezd1=y[3]==-3;kezd2=y[2]==-3;

DSolve[{egy,kezd1},y[x],x]

ds=DSolve[{egy,kezd2},y,x]

y'[3]/.ds

ClearAll[ds];
ds[x_,c_]:=Simplify[DSolve[{egy,y[1]==c},y[x],x][[1,1,2]]]

ds[x,c]

Plot3D[Evaluate[ds[x,c]],{x,0,3},{c,-1,2},AxesLabel->{x,c,y}]

Paraméterekkel is.

DSolve[{y'[x]==a x+b y[x]},y[x],x]

Szakaszonként definált bemenet

DSolve[{y''[x]-y[x]==Max[x,(x^2)],y[0]==1,y'[0]==2},y,x]

Változó együtthatós másodrendű lineáris egyenletek

Euler-féle is.

DSolve[(x^2)y''[x] +4x y'[x]+7y[x]==0,y[x],x]

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

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

Pép

DSolve[{y''[x]-x y[x]==0,y[0]==1,y[9]==1},y,x]

Célbalövéssel

Manipulate[Plot[Evaluate[y[t]/.First[DSolve[ {y''[x]==-x y[x],y[0]==a,y'[0]==b},y,x]]],{t,0,4},Epilog->{Point[{4,1/2}],Green,Arrow[{{0,a},{1,b+a}}],Red,Point[{0,a}]},ImagePadding->25,PlotRange->3],
{{a,1,TraditionalForm[y[0]]},-3,3},
{{b,0,TraditionalForm[y'[0]]},-3,3}]

Lineáris rendszer fázisképe

DSolve[{x'[t]==-2y[t]+3x[t],y'[t]==15x[t]-y[t],x[0]==a,y[0]==b},{x,y}, t];

ParametricPlot[Table[{x[t], y[t]}/.%/.{a->1/(13+m), b->1/(15+m)}, {m,0,20,7}], {t,-2,2},Evaluated->True]

Dinamikusabban

A=(({{-1.1, 0.9}, {-1.4, 0.3}}));
Manipulate[ParametricPlot[Evaluate[MatrixExp[A t,#]&/@pt],{t,0,10},PlotRange->5],{{pt,{{2,0},{0,1},{-3,0}}},Locator},SaveDefinitions->True]

Egy alkalmazás: Görbék rekonstrukciója

T[s_]:=Array[(Subscript[t, #])[s]&,3];
Ν[s_]:=Array[(Subscript[n, #])[s]&,3];
B[s_]:=Array[(Subscript[b, #])[s]&,3]

X[s_]:=Join[T[s],Ν[s],B[s]]

frenet=Flatten[Thread/@{Array[((Subscript[t, #]))'[s]&,3]==κ[s] Ν[s],
Array[((Subscript[n, #]))'[s]&,3]==-κ[s]T[s]+τ[s]B[s],
Array[((Subscript[b, #]))'[s]&,3]==-τ[s]Ν[s]}]

ini=Flatten[Thread/@{Array[(Subscript[t, #])[0]&,3]==UnitVector[3,1],
Array[(Subscript[n, #])[0]&,3]==UnitVector[3,2],
Array[(Subscript[b, #])[0]&,3]==UnitVector[3,3]}]

ClearAll[ini];
ini=Thread[X[0]==Flatten[IdentityMatrix[3]]]

Column[Union[frenet,ini]/.{κ->((Sqrt[2])/2&),τ->((Sqrt[2])/2&)}]

ds=DSolve[Union[frenet,ini]/.{κ->((Sqrt[2])/2&),τ->((Sqrt[2])/2&)},X[s],s]

ParametricPlot3D[Integrate[T[s]/.ds,{s,0,t}],{t,0,8π}]//Timing

ParametricPlot3D[Evaluate[Integrate[T[s]/.ds,{s,0,t}]],{t,0,8π},BoxRatios->{1,1,4},ViewPoint->{0,1,2},PlotStyle->Directive[{Brown,Thickness[0.02]}]]//Timing

m={{1,2},{3,4}};

ArrayFlatten[{{0,0,m},{m,m,0}}]

Coeff[s_]:=ArrayFlatten[{{0,κ [s]IdentityMatrix[3],0},{-κ[s] IdentityMatrix[3],0,τ[s]IdentityMatrix[3]},{0,-τ[s]IdentityMatrix[3],0}}]

Coeff[s]

KroneckerProduct[{{0,κ[s],0},{-κ[s],0,τ[s]},{0,-τ[s],0}},IdentityMatrix[3]]//MatrixForm

Coeff[s]==KroneckerProduct[{{0,κ[s],0},{-κ[s],0,τ[s]},{0,-τ[s],0}},IdentityMatrix[3]]

Last/@frenet==Coeff[s].X[s]

Innen szebb

T[s_]:=Array[(Subscript[t, #])[s]&,3];
Ν[s_]:=Array[(Subscript[n, #])[s]&,3];
B[s_]:=Array[(Subscript[b, #])[s]&,3]

X[s_]:=Join[T[s],Ν[s],B[s]]

Eh[s_]:=KroneckerProduct[{{0,κ[s],0},{-κ[s],0,τ[s]},{0,-τ[s],0}},IdentityMatrix[3]]

Egy=Thread[D[#,s]&/@X[s]==Eh[s].X[s]]

ClearAll[ini];
ini=Thread[X[0]==Flatten[IdentityMatrix[3]]]

ds=DSolve[Union[Egy,ini]/.{κ->((Sqrt[2])/2&),τ->((Sqrt[2])/2&)},X[s],s]

ParametricPlot3D[Evaluate[Integrate[T[s]/.ds,{s,0,t}]],{t,0,8π},BoxRatios->{1,1,4},ViewPoint->{0,1,2},PlotStyle->Directive[{Brown,Thickness[0.02]}]]//Timing

Ez az igazi

ClearAll[GorbeRekonstrualo];
GorbeRekonstrualo[κval_,τval_]:=Module[{X,eh,egy,ini,ds},
X[s_]:=Array[(Subscript[x, #])[s]&,9];
eh[s_]:=KroneckerProduct[{{0,κval,0},{-κval,0,τval},{0,-τval,0}},IdentityMatrix[3]];
egy:=Thread[D[#,s]&/@X[s]==eh[s].X[s]];
ini:=Thread[X[0]==Flatten[IdentityMatrix[3]]];
ds:=DSolve[Union[egy,ini],X[s],s];
ParametricPlot3D[Evaluate[Integrate[Array[(Subscript[x, #])[s]&,3]/.ds,{s,0,t}]],{t,0,8π},BoxRatios->{1,1,4},ViewPoint->{0,1,2},PlotStyle->Directive[{Brown,Thickness[0.02]}]]
]

Cross[{cval,0,gval},#]&/@IdentityMatrix[3]

GorbeRekonstrualo[(Sqrt[2])/2,(Sqrt[2])/2]//MatrixForm

GorbeRekonstrualo[0,1]//MatrixForm

GorbeRekonstrualo[1,0]//MatrixForm

Nemlineáris egyenletek és rendszerek

Egzakt

Ezt akarjuk megoldani:

ex=(4(x^3)((y[x])^3)-2x y[x])+(3(x^4)((y[x])^2)-(x^2))y'[x]==0;

P[x_,y_]:=ex[[1,{1,2}]]/.{y[x]->y}
Q[x_,y_]:=(ex[[1,3]]/.{y[x]->y})/y'[x]

{P[x,y],Q[x,y]}

Keresztbevett parciálisok (+tartomány!!)

D[P[x,y],y]==D[Q[x,y],x]

cpP=ContourPlot[P[x,y]==0,{x,-4,4},{y,-1,1},ContourStyle->Directive[Red,Thick]]

Q[x,y]

cpQ=ContourPlot[-(x^2)+3 (x^4) (y^2)==0,{x,-4,4},{y,-1,1}]

Show[cpP,cpQ](*Tooltip!*)

A primitív függvény meghatározása

A megoldás, F[x,y]=állandó, szintvonalak

F[x_,y_]=Integrate[Q[x,y],y]+Integrate[P[x,y]-Integrate[D[Q[x,y], x], y],x]

Az (1,1) ponton átmenő megoldások:

Union[Solve[F[x,y]==F[1,1],y]]

Plot3D[Evaluate[F[x,y]],{x,-4,4},{y,-1,1}]

(*néhány szintvonal*)
ContourPlot[F[x,y],{x,-5,5},{y,-5,5},Contours->20,ContourShading->False]

ContourPlot[F[x,y],{x,-1.5`,1.5`},{y,-1.5`,1.5`},Contours->20,ContourShading->False,PlotPoints->50]

ContourPlot[F[x,y],{x,-1.5`,0},{y,0,1.5`},Contours->20,ContourShading->False,PlotPoints->50]

Lehetne még: Integráló tényező

Riccati

ric=DSolve[{y'[x]==(x^2)+((y[x])^2)},y[x],x][[1,1,2]];

Together[FunctionExpand[ric]]

Az elemi függvény fogalmának csak történeti jelentősége van. (HF: Program arra, h elemi-e egy fv?)

Nemlineáris rendszer

DSolve[ {y'[x]==Exp[z[x]]+1,z'[x]==y[x]-x},{y,z},x]//Quiet

Implicit egyenlet

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

DAE (Differenciál-algebrai egyenlet)

DSolve[{2 y'[x]+z'[x]==4y[x]+x,y[x]-z[x]==1},{y[x],z[x]},x]

Tyihonov, szinguláris perturbáció, kémiában XX. század legelején, fizikában?

Mire jó még a szimbolikus eredmény?

Adott a gradiens

DSolve[{D[f[x, y], x]==2x (y^3)+y Cos[x y], 
                D[f[x,y],y]==3(x^2)(y^2)+x Cos[x y]},f[x,y],{x,y}]

Adott a képzetes rész

Bizonyítsák be!

f[z_]:=2I v[z/2,-I z/2]

Például

ComplexExpand[(((x+I y))^2)]

v[x_,y_]:=2x y

2I v[z/2,-I z/2]

Generátorfüggvény

GeneratorForStirlingS1[x_,y_] = FullSimplify[u[x, y] /.DSolve[{(1-x)D[u[x,y],x]==y u[x, y],u[0, y]==1},u, {x,y}][[1]], 
   x≥0 ∧ y≥0]

SeriesCoefficient[Series[GeneratorForStirlingS1[x,y],{x,0,10},{y,0,4}],{10,4}]10!

Impulzusválasz

DSolve[y'''[x]-5y''[x]+9y'[x]-5y[x]==DiracDelta''[x]+2DiracDelta'[x]+DiracDelta[x]&&y[-1]==0&&y'[-1]==0&&y''[-1]==0,y[x],x]//FullSimplify

Görbesereg differenciálegyenlete

Eliminate[{y[x]==C (x^3),D[y[x]==C (x^3),x]},{C}]

Eliminate[{(((x-a))^2)+b ((y[x])^2)==1,D[(((x-a))^2)+b ((y[x])^2)==1,x],D[(((x-a))^2)+b ((y[x])^2)==1,{x,2}]},{a,b}]

curvefamilyODE[indep_,eq_,params_]:=Module[{},Eliminate[FoldList[D[#1,indep]&,eq,params],params]]

Rengeteg sok feladat

Philippov 17.

curvefamilyODE[x,y[x]==C (E^(k x)),{C,k}]

Philippov 18. OK, better than in the problem book.

curvefamilyODE[x,y[x]==(((x-C))^3),{C}]

Philippov 19. OK

curvefamilyODE[x,y[x]==C (x^3),{C}]

Philippov 20. OK

curvefamilyODE[x,y[x]==Sin[x+C],{C}]

Eliminate[{y[x]==Sin[x+C],D[y[x]==Sin[x+C],x]},C]

Reduce[{y[x]==Sin[x+C],D[y[x]==Sin[x+C],x]},{C}]

Eliminate[%,C]

Philippov 21. OK

curvefamilyODE[x,(x^2)+C ((y[x])^2)==2y[x],{C}]

Philippov 22. OK

curvefamilyODE[x,((y[x])^2)+C x==(x^3),{C}]

Philippov 23. OK

curvefamilyODE[x,y[x]==C(((x-C))^2),{C}]

Szebb így?

Map[Factor,%,1]

Philippov 24.

curvefamilyODE[x,C y[x]==Sin[C x],{C}]

Upon taking the derivative and the square of the original equation one obtains

(((C y'[x]))^2)==(C^2) ((Cos[C x])^2)

thus

(C^2)==(1-((y'[x])^2))/((y[x])^2.)

Substituting C into the origianl equation one gets

y'[x]==Cos[((x(Sqrt[1-((y'[x])^2)]))/(y[x]))]

How to obtain this automatically?

Philippov 25. OK

curvefamilyODE[x,y[x]==a (x^2)+b (E^x),{a,b}]

Simplify[%]

Hogyan lehet egyenletet egy oldalra rendezni?

Simplify[Subtract@@%]==0

Philippov 26. OK

curvefamilyODE[x,(((x-a))^2)+b ((y[x])^2)==1,{a,b}]

Kijon osztas nelkul is kezzel ?? ??

Philippov 27.

curvefamilyODE[x,Log[y[x]]==a x+b y[x],{a,b}]

FullSimplify[%[[2]]]

y[x]>0 is missing.

Philippov 28. OK

curvefamilyODE[x,y[x]==a (x^3)+b (x^2)+c x,{a,b,c}]

Philippov 29. OK

curvefamilyODE[x,x==a ((y[x])^2)+b y[x]+c,{a,b,c}]

KKM, p. 60, Ex. 1 OK

curvefamilyODE[x, ((x^2)/(a^2))-(((y[x])^2)/1)==1,{a}]

KKM, p. 60, Ex. 2

curvefamilyODE[x, y[x]==a(1-(E^(-(x/a)))),{a}]

a≠0, although we do not construct the whole theory. y'[x]>0

KKM, p. 60, Ex. 3

kkm3=curvefamilyODE[x, x Cos[α]+y[x] Sin[α]-1==0,{α}]

One may, if (s)he wishes, ask for simplifications with not too impressing results.

Simplify[kkm3[[1]]]

kkm3[[2]]

How to obtain this:

sol=Solve[Cos[α]+y'[x] Sin[α]==0,y'[x]]

mmr=Map[Minus,Reverse[sol[[1,1]]]]

Reverse[sol[[1]]]

szab={Cos[α]->((Cot[α])/(Sqrt[1+((Cot[α])^2)])),Sin[α]->(1/(Sqrt[1+((Cot[α])^2)]))}

x Cos[α]+y[x] Sin[α]-1==0/.szab/.Map[Minus,Reverse[sol[[1,1]]]]

FullSimplify[%]

Ezt szeretném:

x y'[x]+y[x]==(Sqrt[1+((y'[x])^(2 ))])

KKM, 238. OK

curvefamilyODE[x, y[x]==a/x,{a}]

KKM, 239. OK

curvefamilyODE[x, (x^2)-((y[x])^2)==a x,{a}]

KKM, 240.

curvefamilyODE[x, y[x]==a (E^(x/a)),{a}]

a≠0, y'[x]>0

KKM, 241.

curvefamilyODE[x, y[x]==C x-C-C (x^2),{C}]

A konyvben:

((y'[x])^2)+y'[x]-x y'[x]+y[x]==0

Egyszerubb? Egyenerteku?

KKM, 242. OK

curvefamilyODE[x, y[x]==(E^x)(a x+b),{a,b}]

Irjuk igy:

d[y_]:=D[y,x]

d[y[x]]

(((d-1))^2)[y[x]]==0

Csak ez mukodjon, de hogyan megy?

KKM, 243.

curvefamilyODE[x, ((y[x])^2)==2C x+(C^2),{C}]

In problem books the equation is divided through by y[x] , excluding some solutions.

KKM, 244. OK

curvefamilyODE[x, y[x]==a (x^2)+b x+c,{a,b,c}]

KKM, 245. OK

curvefamilyODE[x, y[x]==(Subscript[C, 1]) x+((Subscript[C, 2])/x)+(Subscript[C, 3]),{(Subscript[C, 1]),(Subscript[C, 2]),(Subscript[C, 3])}]

KKM, 246.

curvefamilyODE[x, (((x-a))^2)+(((y[x]-b))^2)==1,{a,b}]

If you see a simplifying possibility, you may request it.

Map[Factor,%,1]

KKM, 247. OK

curvefamilyODE[x, y[x]==(Subscript[C, 1])(E^x)+(Subscript[C, 2])(E^(-x)),{(Subscript[C, 1]),(Subscript[C, 2])}]

KKM, 248.

curvefamilyODE[x, y[x]==a Sin[x+α],{a,α}]

We received more possibilities which have nothing to do with the solution; that is given as the first term.

Matveev, 62. OK

curvefamilyODE[x,y[x]==(1/3)(x^3)+C,{C}]

Matveev, 63.

curvefamilyODE[x,y[x]==(Sqrt[1-(x^2)])+C,{C}]

One would be happier to see a warning that x≠1, or that -1<x<1.

Matveev, 64.

curvefamilyODE[x,y[x]==Log[x]+C,{C}]

The warning x>1 would be more appropriate from the mathematical point of view—if one only considers real numbers.

Néhány alapfogalom

Megoldás, trajektória

Inga

nds=NDSolve[{x'[t]==y[t],y'[t]==-Sin[x[t]],x[0]==0,y[0]==1},{x,y},{t,0,15}][[1]];
ParametricPlot3D[Evaluate[{{0,x[t],y[t]},{t,x[t],-1.5},{t,1.5,y[t]},{t,x[t],y[t]}}/.nds],{t,0,15},PlotStyle->{Directive[Thick,Blue],Directive[Thick,Green],Directive[Thick,Brown],Directive[Thick,Red]},BoxRatios->{4,1,1},AxesLabel->{"t","x(t)","y(t)"},PlotLabel->"Az inga egyenletének megoldása,\na megoldás koordinátafüggvényei  és\n trajektóriája a fázissíkon"]

Első integrál

φ[p_,q_]:=((q^2)/2)-Cos[p];f[p_,q_]:={q,-Sin[p]};

D[φ[p,q],{{p,q}}].f[p,q]

ContourPlot[((y^2)/2) - Cos[x], {x, -5, 8}, 
       {y, -2, 2},AspectRatio->Automatic]

Brüsszelátor

brusszel={x'[t]==1+a ((x[t])^2)y[t]-(5+1)x[t],y'[t]==a(5 x[t]-a ((x[t])^2)y[t]),x[0]==8,y[0]==3}/.{a->1.96};

Framed[TraditionalForm[Column[brusszel]]]

ClearAll[nds];
nds=NDSolve[brusszel,{x,y},{t,0,100}][[1]];
ParametricPlot3D[Evaluate[{{90,x[t],y[t]},{t,x[t],0.5},{t,1.6,y[t]},{t,x[t],y[t]}}/.nds],{t,90,100},
BoxRatios->{1.5,1,1},
PlotRange->{{90,100},{0.4,1.6},{0.5,4.}},
AxesLabel->{"t","x[t]","y[t]"},  
PlotStyle->{Directive[Thick,Blue],Directive[Thick,Brown],Directive[Thick,Green],Directive[Thick,Red]},PlotLabel->Style["A Brüsszelátor megoldása,\na megoldás koordinátafüggvényei és\n trajektóriája a fázissíkon","Subsection"],LabelStyle->Directive[Orange,FontFamily->"Helvetica"]]

Iránymező, irányvonalak (elemi kvalitatív vizsgálat)

Ezt vizsgáljuk: y'(x)=2x-y(x), 
azaz (OverTilde[x])'(t)=1     (OverTilde[y])'(t)=2(OverTilde[x])(t)-(OverTilde[y])(t)

sp=StreamPlot[{1,2 x-y},{x,-2,4},{y,-2,4}]

l[k_]:=2 x-k;(*irányvonal, amit k meredekséggel metszenek a megoldások*)
iranyok=Plot[{l[-3],l[0],l[1],l[2]},{x,-3,4},PlotRange->{-3,4},PlotStyle->{{Thickness[0.008`],Black},{Thickness[0.008`],Red},{Thickness[0.008`],Black},{Thickness[0.008`],Black}}]

ClearAll[ds];
ds[x_,c_]:=Simplify[DSolve[{y'[x]==2x-y[x],y[1]==c},y[x],x][[1,1,2]]]

pl[c_]:=Plot[Evaluate[ds[x,c]],{x,-2,4},PlotStyle->{Thickness[0.008],((1+Sign[c])/2)Green+((1-Sign[c])/2)Blue},PlotRange->{-2,4}];

Show[{pl[2],pl[-3]}]

Show[sp,iranyok,pl[1],pl[2],pl[3],pl[-1],pl[-2],pl[-4],PlotRange->{-2,4},Frame->True,Axes->True]

Közelítő módszerek

A pontos megoldás

ClearAll[pontos,y,x];
pontos=DSolve[{y'[x]==x y[x], y[0]==2},y[x],x][[1]]

Euler-módszer (vö. általános tételek)

NestList[ff,x0,2]

euler[f_,h_,x0_,y0_,n_:3]:=Module[{lep,x,y},lep[{x_,y_}]:={x+h,y+h f[x,y]};NestList[lep,{x0,y0},n]]

ClearAll[F];
F[x_,y_]:=x y;

pont=Plot[y[x]/.pontos,{x,0,1},PlotLabel->"pontos"]

abra[n_]:=Module[{},lp=ListPlot[euler[F,((1.`)/n),0,1,n],Joined->True,PlotStyle->{RGBColor[1,0,0]}];pl=Plot[2(E^((x^2)/2)),{x,0,1}];
Show[pl,lp,PlotLabel->"Lépésszám: "<>ToString[n]]]

GraphicsGrid[{{abra[1],abra[100]},{pont}},ImageSize->300]

Szukcesszív approximáció (vö. általános tételek)

Az eszköz

NestList[φ,(Subscript[x, 0]),3]

nl=NestList[Cos,-6.7,12]

ListPlot[nl,PlotStyle->{PointSize[0.02`]}]

Alkalmazzuk a differenciálegyenletre

Az egyenértékű integrálegyenlet:

(∫_0^x)y'[ξ]\[DifferentialD]ξ==(∫_0^x)ξ y[ξ]\[DifferentialD]ξ/.{y[0]->1}

Az operátor, amit iterálunk:

ClearAll[A];
A[φ_]:=1+Integrate[x φ, x]

nl=NestList[A , 1,12]

Összehasonlítás

ClearAll[pontos,y,x];
pontos=DSolve[{y'[x]==x y[x], y[0]==2},y[x],x][[1]]

Plot[Evaluate[{y[x]/.pontos,nl[[3]]}],{x,0,4},PlotStyle->{Red,Blue}]

Találjuk ki a hatványsor összegét

Becslések az értelmezési tartományra?

És ha NDSolve-ot használunk?

Rendszerekre?

Maradék

DynamicModule[{v1={2,0},v2={-1,1}},Dynamic@Graphics[{Circle[],Red,GeometricTransformation[Circle[],Transpose[{v1,v2}]],Green,Line[{{0,0},v1}],Line[{{0,0},v2}],Locator[Dynamic[v1]],Locator[Dynamic[v2]]},PlotRange->3]]

Házi feladatok

1. Az előző órán tanult becslési módszerek bármelyikét alkalmazzák olyan adatokra, amelyek 1, vagy saját mérésből 2, vagy a hálóról 3. vagy a Mathematica ?*Data adataiból származik.

2. Keressenek a hálón Mathematica egysorosokat (oneliner), és próbálják megérteni a működésüket.

3. Mit tudnak közönséges differenciálegyenletek periodikus megoldásainak létezéséről? Ez lesz a téma legközelebb.