MATH529 Lesson02

First-order differential equations

Exact equations

The equation z(x,y)=c is an implicit definition of a curve C and the differential dz= M dx+N dy  evaluates to zero on the curve  (M=(∂_x)z, N=(∂_y)z). This implies that  =c is a solution to the ODE
((dy)/(dx))=-(M/N), or y'(x)=-N/M=f(x,y) . A solution to the ODE exists if  is continuous and differentiable, hence (x,y) has continuous second derivatives, in which case (Subscript[z, x y])=(Subscript[z, y x]).

Criterion for an exact differential

In general, the differential form M(x,y)dx+N(x,y) dy is said to be exact if there exists some z(x,y) such that (Subscript[z, x])=M(x,y), (Subscript[z, y])=N(x,y). When  is twice differentiable, the condition for an exact differential form is 
(Subscript[M, y])=(Subscript[N, x]) or (∂_y)M=(∂_x)N.

Example 1: Exact differential corresponds to ODE solution

Curves defined by z(x,y)=c .

z[x_,y_] = (x^3)-5x y-(y^2)

cplt=ContourPlot[z[x,y],{x,0,2},{y,-1,5},Contours->Table[c,{c,0,2,0.2}],ContourShading->None]

y'(x)=-N/M=f(x,y)

Construct the ODE y'(x)=-N/M=f(x,y), choose an initial condition and superimpose on the above contours

zM[x_,y_]=(∂_x)z[x,y]

zN[x_,y_]=(∂_y)z[x,y]

y'(x)=-N/M=f(x,y)

f[x_,y_]=-zM[x,y]/zN[x,y]

sol=DSolve[{y'[x]==f[x,y[x]],y[0.5]==-1},y[x],x][[1,1]]

splt=Plot[y[x] /. sol,{x,0.5,2},PlotStyle->Red]

Show[{cplt,splt}]

Example 2: Verification of exact differential form

Consider the differential form ((e^(2y))-y cos x y)dx+(2x (e^(2y))-x cos x y+2 y)dy

zM[x_,y_]=(E^(2y))-y Cos[x y]; zN[x_,y_]=2x (E^(2y))-x Cos[x y]+2 y;

Verify that it is exact

(∂_y)zM[x,y]==(∂_x)zN[x,y]

Solving exact differential equations

Given  M(x,y)dx+N(x,y) dy=0, the function z(x,y) that satisfies (Subscript[z, x])=M(x,y), (Subscript[z, y])=N(x,y), is recovered by:
1. Integration with respect to x holding y constant
z(x,y)=∫M(x,y)dx+g(y)

2. To find g(y), differentiate result with respect to y, and set equal to N(x,y)
(∂_y)z=N(x,y)=(∂/(∂y))∫M(x,y)dx+g'(y)

Integrating the result gives g(y)
g(y)=∫N(x,y)dy-∫M(x,y)dx

Example 3: Solving an exact differential equation

Consider the differential equation ((e^(2y))-y cos x y)dx+(2x (e^(2y))-x cos x y+2 y)dy=0

zM[x_,y_]=(E^(2y))-y Cos[x y]; zN[x_,y_]=2x (E^(2y))-x Cos[x y]+2 y;

It is exact

(∂_y)zM[x,y]==(∂_x)zN[x,y]

1. Integration with respect to x holding y constant
z(x,y)=∫M(x,y)dx+g(y)

z[x_,y_]=Integrate[zM[x,y],x]+g[y]

2. Find g(y)=∫N(x,y)dy-∫M(x,y)dx

g[y_]=Integrate[zN[x,y],y]-Integrate[zM[x,y],x]

Since g(y) is now defined, it will be substituted in the definition of z(x,y)

z[x,y]

Verify

{D[z[x,y],x]==zM[x,y],D[z[x,y],y]==zN[x,y]}

Solution by substitutions

In y'=f(x,y) the dependent variable y(x) may be replaced by y=g(x,u).

Homogeneous functions

A function f(x,y) is homogeneous if f(s x, s y)=(s^n) f(x,y) with n the degree of the homogeneous function. A differential equation  M(x,y)dx+N(x,y) dy=0 where M(x,y), N(x,y) are homogeneous functions of the same degree can be reduced to a separable equation by the substitution y= u x.
 M(x,y)dx+N(x,y) dy=M(x,u x)dx+N(x,u x) d(u x)=(x^n)[M(1,u)dx+N(1,u) d(u x)]=0⇒ 
 [M(1,u)+u N(1,u)]dx+x N(1,u) du=0⇒
((dx)/x)+((N(1,u) du)/(M(1,u)+u N(1,u)))=((dx)/x)+F(u) du=0⇒ln x+∫F(u)du=ln c⇒x=c exp[-∫F(u)du]=G(u)

Example

Solve ((x^2)+(y^2)) dx+((x^2)-x y)dy=0

zM[x_,y_]=(x^2)+(y^2); zN[x_,y_]=(x^2)- x y;

F[u_]=Simplify[((zN[1,u])/(zM[1,u]+u zN[1,u]))]

Integrate[F[u],u]

G[u_]=Exp[Integrate[-F[u],u]]

sol = x== c G[y/x]

Bernoulli’s equation

In y'+p(x)y=f(x)(y^n), replace u=(y^(1-n)) to obtain a linear ODE

Example

Solve x y'+y=(x^2) (y^2)

DE = x y'[x] + y[x] == (x^2) (((y[x]))^2)

Substitute u=(y^(1-n))=(y^(1-2))=(y^(-1))

sub = y[x]-> 1/u[x]

Differentiate on both sides of substitution

dsub = Simplify[D[sub,x]]

Use substitution in original equation

uDE =Assuming[{n>1,u[x]>0}, FullSimplify[DE /.{sub,dsub}]]

The DE is now linear

usol=DSolve[uDE,u[x],x][[1,1]]

Verify original equation

y[x_] = 1/u[x] /. usol

Simplify[DE]

Linear Models

In many models the rate of growth is proportional to the current state

Bacterial growth

Let  denote current population,  the growth rate

DE = P'[t] - k P[t] == 0; IC = P[0] == P0;
sol = DSolve[{DE,IC},P[t],t][[1,1]]

Define some common plot options for reuse

pltopt = {GridLines->Automatic, PlotLabel->title, Frame->True, FrameLabel->xylabel};

Plot the solution

Plot[P[t] /. sol /. {k->0.2,P0->1},{t,0,10}, Evaluate[pltopt /. {title->"Population",xylabel->{"t", "P"}}]]

Radioactive decay

Let  denote current number of isotopes,  the decay rate

DE = A'[t] + k A[t] == 0; IC = A[0] == A0;
sol = DSolve[{DE,IC},A[t],t][[1,1]]

Plot[A[t] /. sol /. {k->0.2,A0->(10^10)},{t,0,10}, Evaluate[pltopt /. {title->"Nr. isotopes",xylabel->{"t", "A"}}]]

LR circuit

An in-series electrical circuit of source , resistor , and inductor  gives DE

Plot[A[t] /. sol /. {k->0.2,A0->(10^10)},{t,0,10}, Evaluate[pltopt /. {title->"Nr. isotopes",xylabel->{"t", "A"}}]]

DE = L i'[t] + R i[t] == e[t]; IC = i[0]==0;
sol = DSolve[{DE,IC} /. e[t] -> 110 Cos[100π t],i[t],t][[1,1]]

Plot[i[t] /. sol /. {R->(10^5),L->500},{t,0,0.1}, Evaluate[pltopt /. {title->"RL circuit current",xylabel->{"t", "i"}}],PlotRange->All]

After initial transient, the circuit stabilizes at current imposed by source.

Nonlinear models

Variable population growth rate (logistic equation)

Environmental constraints on population growth can be modeled by the logistic equation



where  is the growth rate when , and  is maximum population

DE = P'[t] == r (1-P[t]/K) P[t]; IC = P[0]==P0;
sol=DSolve[{DE,IC},P[t],t][[1,1]]

Growth from population less than environment carrying capacity

Plot[P[t] /. sol /. {K->(10^5),r->2,P0->100},{t,0,10}, Evaluate[pltopt /. {title->"Logistic population model",xylabel->{"t", "P"}}],PlotRange->All]

Decrease when population is greater than

Plot[P[t] /. sol /. {K->(10^5),r->2,P0->(10^6)},{t,0,10}, Evaluate[pltopt /. {title->"Logistic population model",xylabel->{"t", "P"}}],PlotRange->All]

Multiple initial conditions rendered on same plot

Plot[Table[P[t] /. sol /. {K->(10^5),r->2,P0->(10^p)},{p,1,6}],{t,0,10}, Evaluate[pltopt /. {title->"Logistic population model",xylabel->{"t", "P"}}],PlotRange->All]

Systems of first-order differential equations

Predator-prey model

Foxes  prey on rabbits , but lacking prey the fox population would decrease. The rabbit population increases (vegetation is freely available) in absence of foxes. The maximum number of rabbit- fox encounters is . The resulting (Lotka-Volterra) model is

DE = {f'[t]==-a f[t] + b f[t] r[t], r'[t] == c r[t] - d f[t] r[t]};
IC = {f[0]==4, r[0]==4};
IVP = Flatten[{DE,IC}];
params = {a->0.16, b ->0.08, c -> 4.5, d -> 0.9};

An analytical solution of the above non-linear system is tedious to obtain. Compute a numerical approximation

sol=NDSolve[IVP /. params,{f[t],r[t]},{t,0,20}][[1]]

Plot[{f[t] /. sol,r[t] /. sol},{t,0,20},Evaluate[pltopt /. {title->"Lotka-Volterra model",xylabel->{"t", "f,r"}}],PlotLegends->{"Foxes","Rabbits"}]