Solving Differential Equations

An ordinary differential equation (ODE) is an equation involving a function y(x) and its derivatives y’(x), y’’(x), ... , (y^((n)))(x). The order of the ODE is the highest order of the derivatives occurring. For example, the ODE y’’’(x)+2(e^x)y’’(x)+y(x)y’(x)=(x^4) has order 3. A function y=y(x) which verifies this differential equation is called a solution of the ODE, or an integral curve.

Example 1:  
Solve the ODE  (dy/dx)=((x^2)/(1-(y^2))) and sketch its solutions.

This ODE is of a special type and can be solved using separation of variables. We first write:
(1-(y^2))dy = (x^2)dx 
and integrate each side with respect to y and x.
∫(1-(y^2))dy = ∫(x^2)dx 

The general solution is found (by integration) to be  
y-((y^3)/3)+c=((x^3)/3), where c is a constant parameter.

The implicit curves
 y-((y^3)/3)+c=((x^3)/3) are called the integral curves of the differential equation.

Example 2 (Initial value problem):
Solve the ODE  (dy/dx)=((x^2)/(1-(y^2))), with initial condition y(0)=-1.1

Otherwise said, we need to find the unique solution of the ODE that passes through the point (x0,y0)=(0,-1.1)

Given a point (x0,y0), to find the solution that passes through the point (x0,y0) means to solve the equation (Subscript[y, 0])-(((Subscript[y, 0])^3)/3)+c=(((Subscript[x, 0])^3)/3) and find the value of the parameter c.

Solve[y0=-1.1;x0=0; y0-y0^3/3+c==x0^3/3, c]

Vector field plot

Pick any point P=((Subscript[x, 0]), (Subscript[y, 0])) and compute the tangent line at P to the solution curve passing through the point P. The slope of the tangent line at P is (Subscript[m=(dy/dx)|, P])=(((Subscript[x, 0])^2)/(1-((Subscript[y, 0])^2))), so the equation of the tangent line at P is y-(Subscript[y, 0]) = m(x-(Subscript[x, 0])).

VectorPlot[{vx, vy}, {x, xmin, xmax}, {y, ymin, ymax}]
generates a vector plot of the vector field as a function of x and y.

VectorPlot[{1,x^2/(1-y^2)},{x,-3,3},{y,-3,3},VectorScale->{Tiny,Small, None}]

The option VectorScale determines the length and arrowhead size of field vectors that are drawn.

VectorPlot[{1,x^2/(1-y^2)},{x,-3,3},{y,-3,3},VectorScale->{Tiny,Small, None},StreamPoints->Coarse,StreamColorFunction->Hue]

Stream Plot

StreamPlot[{vx, vy}, {x, xmin, xmax}, {y, ymin, ymax}]
generates a stream plot of the vector field as a function of x and y.

StreamPlot[{1,x^2/(1-y^2)},{x,-3,3},{y,-3,3}]

Plot the integral curves of the ODE (dy/dx)=((x^2)/(1-(y^2))), with the integral curve that passes through the point (0, -1.1) plotted in red.

StreamPlot[
{1,x^2/(1-y^2)},{x,-3,3},{y,-3,3},
StreamScale->Tiny,
StreamPoints->{{{{0,-1.1},Red}, Automatic}}]

Exact Solutions

DSolve[eqn,y,{x,xmin,xmax}]
solves a differential equation for the function y, with the independent variable x between xmin and xmax.

DSolve returns results as lists of rules. This makes it possible to return multiple solutions to an equation. For a system of equations, possibly multiple solution sets are grouped together inside curly brackets { }. You can use the rules to substitute the solutions into other calculations.

It may not always be possible to compute exact solutions of an ODE.

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

Numerical Solutions

NDSolve[eqns, y, {x, xmin, xmax}]
finds a numerical solution to the ordinary differential equations eqns for the function y with the independent variable x in the range xmin to xmax. The output of NDSolve is a list of transformation rules.

sol1=NDSolve[
{y'[x]==x^2/(1-y[x]^2),y[0]==-1.1},y,{x,0,5}]

The solution to the initial value problem from Example 2 is unique. In this case, the output of NDSolve is a nested list with one transformation rule. The transformation rule describes y as a (numerically computed) function of x. Conceptually, think of the transformation rule as something like y → f[x_ ]

Recall : To apply a transformation rule (or a set of transformation rules) to an expression, we use the Replacement operator “/.” (slash-dot).

y[2] /. sol1

We can also plot the function returned by NDSolve[..], just as we plot any other function in the Wolfram language, using the command Plot[..]. Compare the graph below to the stream line that we had plotted in red in the last StreamPlot[..].

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

The command Flatten[..] removes one set of curly brackets.

Flatten[sol1]
y[2] /. Flatten[sol1]

If we plan to use the function returned by NDSolve several times, it makes sense to give it a name, so we do not have to apply the replacement operator "/." (slash - dot) every time.

f[z_]:= Module[{sol0,x,y},
sol0= NDSolve[{y'[x]==x^2/(1-y[x]^2),y[0]==-1.1},y,{x,0,5}];
y[z] /. Flatten[sol0]
]

f[2]

Plot[f[z],{z,0,3},AxesLabel->{"x-axis","y-axis"}, PlotStyle->{Red, Thick}]

Exercise: Study the ODE  (dy/dx)=((4x-(x^3))/(4+(y^3)))

Plot the vector field using StreamPlot[...] in the region x ∈ [-3.5, 3.5], y ∈ [-3.3, 3.3]. On the same plot, plot the trajectory starting at the point (0, 1) in Red, and the trajectory starting at the point (1, -1.5) in Green.

Find the discontinuities of the vector field (where the direction of the arrows changes, and the slopes of the integral curves are ∞ or -∞) using NSolve[expr,vars]. This command attempts to find numerical approximations to the solutions of the system expr of equations or inequalities for the variables vars.

Use NDSolve[...] to solve the ODE from example 2, (dy/dx)=((4x-(x^3))/(4+(y^3))), with initial condition y(0)=1. The range for the variable x is assumed to be [-3,3]. Plot the function that NDSolve returns, on the domain [-3,3].

This ODE is a separable differential equation, that we can solve and get the family of solutions defined by the implicit equation  4y+(1/4)(y^4)=2(x^2)-(1/4)(x^4)+c . Plot together on the same picture the implicit solutions corresponding to the following choices of the constant c=0,1,2,-1,-4,-5, using the function ContourPlot[...]. Is the picture consistent to the Vector Field pictures obtained in part 1 and part 2?