Mathematica  és differenciálegyenletek 2.

Legyenek meg itt is az új házi feladatok.

Házi feladatok

Tanulmányozzuk a logisztikus egyenletet azzal a feltétellel, hogy szüretelünk is, vagy másképp lehalásszuk a halak egy részét, mégpedig állandó sebességgel, azaz tekintsük ezt az egyenletet:

x'[t]==x[t](1-x[t])-h

ahol h pozitív állandó. Ennek értékétől hogyan függ a rendszer viselkedése?

Mi a helyzet periodikus lehalászás esetén, azaz itt:

x'[t]==5x[t](1-x[t])-0.8(1+Sin[2π t])

Rajzoljuk meg a fázisegyenest az alábbi egyenlethez:

x'[t]==x[t](1-x[t])(x[t]-4)

Rajzoljuk meg a fázisegyenest az alábbi egyenlethez:

α'[t]==A Exp[-(E/(R ((Subscript[T, 0])+β t)))](((1-α[t]))^n)

Itt minden együttható pozitív, kivéve n -et, ez lehet 0,1,2,3.

Rajzoljuk meg a fázissíkot kétváltozós lineáris rendszerek összes kvalitatíve különböző esetére. Tehát az egyenlet:

{x'[t]==a x[t]+b y[t],y'[t]==c x[t]+d y[t]}

és az A={{a,b},{c,d}} mátrix sajátértékei valósak, különbözők, egybeesők, komplexek stb.

Kérdések és megjegyzések a Mathematica Workoshop után

Legelőször ismertettem a Wolframéknak küldött levelemet; idemásolom.

Questions and remarks

Is it possible to Import data into a notebook using CDF? Or any other data fro outside?

What is the real difference between the Student license and the full one?

I try to give my lectures in such a way that I am using a client machine and the notebook is reached through the local network. When I use Animate or Manipulate sometimes I cannot stop it. A related problem may be that scrolling can also not be halted.

During your talk you used a magnifier. Is it a part of your operating system or is it a part of Mathematica?

Valaki elárulta (BG?), hogy ilyen van a Windowsban is, tudja m ár valaki rendesen használni?

It takes usually very long time if I would like to acquire data from the Wolfram server using any of the functions *Data.

If I would like to get information about something what is not an explicit function name than I may get no answer at all, even if the word I have typed in is an Option name. Is this the case you mentioned that you do not really care because there are more important things to do then polishing the help?

Import understands files with the extension .net as well. An important fact (missign from the help) nowadays when so many people are crazy about networks and graphs.

Functions

This is how I teach them.

Function[{x}, body]

corresponds exactly to the mathematical form x↦body. Or, if the domain of the function is also given (as it should always be the case): ℝ∋x↦body.

Everyone has learned the identity function of a nonempty set A: (Subscript[id, A])(x):=x  for all x∈A.  (One may write id—without the index—if the underlying set is the set of real or complex numbers.) This corresponds exactly to the Mathematica form

#&[ff]

Similarly, if (Subscript[pr, 1]):X×Y⟶X is defined to be (Subscript[pr, 1])(x,y):=x, then one sees that this function is the same as

#1&[p,k]

in Mathematica. It would make life easier for the users if we would not say nebulous sentences about pure functions, slots, anonymous functions (let alone lambda calculus!) etc. The above functions are taught at the very beginning of calculus (then they are forgotten also within mathematics...).

Mátrix elemeinek módosítása

Azután arról a feladatról volt szó, amelynél me kell keresni egy mátrixban bizonyos elemek helyét, majd egy másik mátrixban az ezen helyeken levő elemek értékét alkalmasan felül kell írni. A megoldás több szempontból is tanulságos.

Első megoldás

Az első mátrixban nézzük meg, hol van a betű.

A={{a,b,e},{e,r,a},{s,a,s}};

pos=Position[A,a]

A másik mátrix, amelynek elemeit módosítani akarjuk:

K=Array[(Subscript[k, #1,#2])&,{3,3}];
MatrixForm[K]

Válasszunk ki egy alkalmas elemet (fejcsere!)

Apply[Part[K,#1,#2]&,pos[[1]]]

Vérszemet kaptunk

Apply[Part[K,#1,#2]&,pos,1]

(A harmadik argumentum azt adja meg, hogy melyik szinten akarjuk végrehajtani a fejcserét.)

Apply[Part[K,#1,#2]&,pos,1]={α,β,γ}

Azt hiszi a program, hogy a (beépített, védett: Protected) Apply függvény definícióját akarjuk felülírni. Ha ezt akarnók, akkor az Unprotect[Apply] utasítást kellene kiadnunk. De nem ez a megoldás.

Evaluate[Apply[Part[K,#1,#2]&,pos,1]]={α,β,γ};
K//MatrixForm

Most alkalmazzuk az előzőket K helyett egy olyan mátrixra, amelyiknek minden eleme 0. Másutt megszoktuk, hogy egy mátrixot föltöltünk nullákkal, majd módosítgatjuk. Mathematicában se megy az, hogy egy létre nem hozott mátrix elemeit töltjük föl:

Evaluate[Apply[Part[KK,#1,#2]&,pos,1]]={α,β,γ};

Szóval nézzük a nulla mátrixot.

Evaluate[Apply[Part[Nulla=ConstantArray[0,{3,3}],#1,#2]&,pos,1]]={α,β,γ}

Nulla

A baj hasonló a fentihez: a nulla számnak szerettük volna értéket adni. De a megoldás nagyon eltér a fentitől. Definiálunk egy ritka mátrixot (ennek minden nem definiált eleme nulla, a nullától különbözők a Rule függvénnyelkapnak értéket, tehát szükségünk van egy sor szabályra.)

Rule[#,1]&/@pos

sa=SparseArray[Rule[#,1]&/@pos,{3,3}]

MatrixForm[sa]

Szeretnénk tudni, hány egyes van a mátrixban.

Total@Total@sa

Count[sa,1,2]

Ez miért nem megy, amikor erre remekül működik?

Count[A,a,2]

Count[Normal@sa,1,2]

Tanulság: Ugyan általában a ritka mátrixokkal a program úgy számol, mint a közönségesekkel:

Eigensystem[sa]

néha azonban (mikor??) a Normal függvény bevetésével nógatni kell. Egy másik alkalmazás (haladóknak!):

gc=GraphicsComplex[{{1.2,3.5},{5.6,7.2}},{Point[1],Line[{1,2}]}];

g=Normal[gc]

{Graphics[gc],Graphics[g]}

Mathematica  és differenciálegyenletek 2.

Az előző házi feladatok megoldásáról Holló Gábor és  Enkhbold Enkhbyamba megoldásai alapján

Amit az óra elején csináltunk, azt reprodukáljuk a ParametricNDSolve(Value) függvényekkel.

DSolve[y[x]+y'[x]+y''[x]==a(Sin[bx]+Cos[bx]),y[x],x]

Tanulság: egymás mellé írt betűk új azonosítót jelentenek.

Amit az óra elején csináltunk, azt reprodukáljuk a ParametricNDSolve(Value) függvényekkel.

ClearAll[egy,kezd,ds];
egy=y'[x]==a x-y[x];
kezd=y[3]==-3;

ds=ParametricNDSolve[{egy,kezd},y,{x,0,10},a];

y/.ds

y1=y[2]/.ds

y1[10]

y[2][10]/.ds

Plot3D[Evaluate[y[a][x]/.ds],{a,-1,2},{x,0,10},AxesLabel->{Style[a,24,Bold],Style[x,24,Bold],Style[y,24,Bold]}]

Plot[Evaluate[Table[y[a][x]/.ds,{a,-1,2,0.1}]],{x,0,10},AxesLabel->{x,y}]

Evaluate[y[a][x]/.ds]

y[2][10]/.ds

Evaluate[y[1.1][3.3]/.ds]

Manipulate[Plot[y[a][x]/.ds,{x,0,10},AxesLabel->{x,y},PlotRange->{-40,40}],{a,-3,2,0.1},SaveDefinitions->True]

Most úgy látom, hogy a ParametricNDSolve eredményét nem leht manipulálni, ez biztosan nem igaz.

Milyen a érték mellett válik monotonná?

ds=DSolve[{egy,kezd},y,x][[1]]//Simplify

y'[x]/.ds//Simplify

Solve[(3+2 a) (E^3)+a (E^x)==0,a]

Ha a≥0, akkor nincs megfelelő x.  Ha a<0, akkor

Solve[(3+2 a) (E^3)+a (E^x)==0,x,Reals]

Hoppá, -3/2<a<0 esetén van a deriváltnak nullahelye.

RegionPlot[(3+2 a) (E^3)+a (E^x)>0,{a,-3,1},{x,0,4},FrameLabel->{Style[a,Bold,Red,24],Style[x,Bold,Blue,24]}]

sym=DSolve[{egy,kezd},y[x],x]

Plot[Evaluate[y[x]/.sym/.a->-3/2],{x,0,10},AxesLabel->{x,y},PlotRange->{-40,40},PlotLabel->Simplify@Evaluate[y[x]/.sym/.a->-3/2]]

Az Escherichia coli baktérium 30 percenként osztódik. Egy példányát, amelynek tömege 2 (10^(-12)) g, beléhelyezzük egy igen nagy kémcsőbe/Petri-csészébe, és ellátjuk a szükséges táplálékkal. Mennyi idő múlva éri el a baktériumpopuláció össztömege a Földét (amely 5.9742 (10^27) g)?

baci=First@DSolve[{m'[t]==k m[t],m[0]==2 (10^(-12))},m,t]

sol=Solve[((m[t+30]/m[t])/.baci)==2//Simplify,k,Reals][[1]]

tsol=First@Solve[(m[(SubStar[t])]/.baci/.sol)==5.9472 (10^27),(SubStar[t])]//Quiet

(((SubStar[t])/.tsol)/1440)

Nem kell három nap se.

Írjuk fel a sorosan kapcsolt ellenállásból, kondenzátorból és tekercsből álló RCL-körben folyó áram időbeli változását leíró másodrendű, lineáris, állandó együtthatós differenciálegyenletet, és tanulmányozzuk a megoldások viselkedését a paraméterek függvényében.

(* Lq''+Rq'+1/Cq=E(t) *)
ClearAll[ds]
ds[t_,L_,R_,CC_]=Simplify[DSolve[{L q''[t]+R q'[t]+1/CC q[t]==Sin[t],q[0]==0,q'[0]==1},q[t],t][[1,1,2]]]

Manipulate[Plot[Evaluate[ds[t,L,R,CC]],{t,0,200},AxesLabel->{t,q},PlotRange->{-10^-5,10^-5},PerformanceGoal->"Quality"],{L,0.1,1,0.01},{R,0.1,10,0.01},{CC,10^-6,10^-5,10^-7}]

Más szavakkal: tekintsünk egy testet, amelyre közegellenállás és rugóerő hat, írjuk föl az időbeli változását leíró másodrendű, lineáris, állandó együtthatós differenciálegyenletet, és tanulmányozzuk a megoldások viselkedését a paraméterek függvényében.

Tanulmányozzuk az E+S<=>C->E+P  Michaelis—Menten-reakció kinetikai differenciálegyenletét: e'=-(Subscript[k, 1])e s+(Subscript[k, -1])c+(Subscript[k, 2])c
s'=-(Subscript[k, 1])e s+(Subscript[k, -1])c
c'=+(Subscript[k, 1])e s-(Subscript[k, -1])c-(Subscript[k, 2])c
p'=+(Subscript[k, 2])c.

A szokásos kezdeti feltételek: e(0)=(Subscript[e, 0])>0   s(0)=(Subscript[s, 0])>>(Subscript[e, 0])     c(0)=0     p(0)=0.

ClearAll[ds]
k1=100;
k1r=1;
k2=0.1;
ds=NDSolve[{e'[t]==-k1 e[t] s[t]+k1r c[t]+k2 c[t],s'[t]==-k1 e[t] s[t]+k1r c[t],c'[t]==k1 e[t] s[t]-k1r c[t]-k2 c[t],p'[t]==k2 c[t],e[0]==0.01,s[0]==1,c[0]==0,p[0]==0},{e,s,c,p},{t,0,1500}]

Plot[Evaluate[{e[t],s[t],c[t],p[t]}/.ds],{t,0,1500},PlotLegends->{e[t],s[t],c[t],p[t]},AxesLabel->{t,c}]

Érdekes, hog majdnem törést mutatnak a görbék. Jó lenne valós paramétereket szerezni.

Plot[Evaluate[{c[t],c'[t]}/.ds],{t,0,1000},PlotLegends->{c[t],c'[t]},AxesLabel->{t,c}]

ClearAll[ds,k1,k1r,k2]
ds=ParametricNDSolve[{e'[t]==-k1 e[t] s[t]+k1r c[t]+k2 c[t],s'[t]==-k1 e[t] s[t]+k1r c[t],c'[t]==k1 e[t] s[t]-k1r c[t]-k2 c[t],p'[t]==k2 c[t],e[0]==0.01,s[0]==1,c[0]==0,p[0]==0},{e,s,c,p},{t,0,1500},{k1,k1r,k2}]

c[1,1,1][10]/.ds

Manipulate[Plot[Evaluate[c[k1,k1r,k2][t]/.ds],{t,0,1000},AxesLabel->{t,c},PlotRange->{0,0.02}],{k1,1,1000,1},{k1r,0.1,100,0.1},{k2,0.01,1,0.01}]

Házi feladatok 2.

Ezeket inkább majd az órán beszéljük meg, közösen, de szabad gondolkodni rajtuk.

Oldjuk meg (papíron, ceruzával) a

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

Euler-egyenletet az ,,x=(e^t)" változó bevezetésével. Mit is jelent ez? ReplacementVariables.m: változócseréhez, de a diszkussziót nem végzi el.

Names["ReplaceVariables`*"]

ReplaceVariables[{7 y[x]+4 x y'[x]+(x^2) y''[x]==0},t->Log[x],y[x]->Y[t],Y[t]]

Írjunk programot a célbalövésre a

{y''[x]==-x y[x],y[0]==2,y[4]==1/2}

feladat numerikus megoldására.

pnds=ParametricNDSolve[{y''[x]==-x y[x],y[0]==2,y'[0]==a},y,{x,0,4},{a}]

ff[a_]:=(y[a][4]/.pnds)

Plot[ff[a],{a,-2,3}]

fr=FindRoot[ff[a]==1/2,{a,1}]

ff[a]/.fr//Chop

ClearAll[celba];celba[egy_,perem1_,perem2_,Tmin_,Tmax_,kezdseb1_:1]:=Module[{},pnds=ParametricNDSolve[{egy,y[0]==perem1,y'[0]==a},y,{x,0,Tmax},{a}];
ff[a_]:=(y[a][Tmax]/.pnds);
Plot[ff[a],{a,-2,3},AxesLabel->(Style[#,Bold,16,FontFamily->"Helvetica"]&/@{"a","ff(a)"}),PlotLabel->First@FindRoot[ff[JóKezdetiSebesség]==perem2,{JóKezdetiSebesség,1}]]]
celba[y''[x]==-x y[x],2,1/2,0,4]

y[0.7][4]/.%

IntegerPart

{y''[x]==-x y[x],y[0]==2,y[4]==1/2}

Írjunk programot az integráló tényező meghatározására.

DSolve[y'[x]==-((2x y[x]+((y[x])^3))/(3(x^2)+5x(((y[x]))^2))),y[x],x]

nex=(2x y[x]+((y[x])^3))+y'[x](3(x^2)+5x(((y[x]))^2))==0

P[p_,q_]:=nex[[1,{1,2}]]/.{y[x]->q,x->p}
Q[p_,q_]:=(nex[[1,3]]/.{y[x]->q,x->p})/y'[p]

{P[p,q],Q[p,q]}

D[P[p,q],q]==D[Q[p,q],p]

D[μ[q]P[p,q],q]==D[μ[q]Q[p,q],p]

D[μ[q]P[p,q],q]-D[μ[q]Q[p,q],p]==0//Simplify

ds=First@DSolve[-2 μ[q]+q μ'[q]==0,μ,q]/.{C[1]->1}

D[μ[q]P[p,q],q]==D[μ[q]Q[p,q],p]/.ds//Simplify

int1[nex_]:=Module[{},P[p_,q_]:=nex[[1,{1,2}]]/.{y[x]->q,x->p};
Q[p_,q_]:=(nex[[1,3]]/.{y[x]->q,x->p})/y'[p];If[D[P[p,q],q]-D[Q[p,q],p]===0,Return["exact"]];
rule=First@Solve[Simplify[D[μ[q]P[p,q],q]-D[μ[q]Q[p,q],p]==0],μ'[q]];ds=DSolve[Equal@@First@rule,μ[q],q]/.{C[1]->1};
{μ[q]P[p,q],μ[q]Q[p,q]}/.ds];

int1[(2x y[x]+(((y[x]))^3))+y'[x](3(x^2)+5x (((y[x]))^2))==0]

Ées már jöhet az egzaktmegoldó ☺

Programozott fűtés

nds=NDSolve[{α'[t]==(E^(-1/(273+t)))(1-α[t]),α[0]==1.4},α[t],{t,0,25}];

Plot[α[t]/.nds,{t,0,25}]