Mathematica notes for Section 9.2

Examine the DE

(dy/dx)=((1+x)/(x y))

Examine the IVP

(dy/dx)=((1+x)/(x y)), y(-3.5)=0.1

Note you can solve the DE by hand since it separable. I am going to use Mathematica to do everything so you can see how Mathematica would work.

If you like, you can choose Evaluation > Evaluate Notebook from the menu to evaluate all the cells.

Direction Fields: Global Picture of Solution to DE

Define the function which appears in our differential equation dy/dx=f(x,y):

f[x_,y_] = (1+x)/(x*y)

Here is how you plot the direction field:

field=VectorPlot[{1,f[x,y]},{x,-4,4},{y,-4,4},VectorScale->{Tiny,Small,None},BaseStyle->16,FrameLabel->{"x","y"}]

You can add streams to the vector field, which are just some members of the family of curves that satisfy the DE.

fieldstreams=VectorPlot[{1,f[x,y]},{x,-4,4},{y,-4,4},VectorScale->{Tiny,Small,None},BaseStyle->16,FrameLabel->{"x","y"},StreamPoints->30,StreamStyle->Purple,StreamScale->Large]

Notice the arrows on the streams change direction at y=0.

Euler's Method: Numerical Solution to IVP

Here is how you implement Euler's method. The semicolon at the end of the line surpresses output.

Clear[x,y]
f[x_,y_]=(1+x)/(x*y);
x[0]=-2.1;
y[0]=0.1;
h = 0.05;
x[n_]:=x[n]=x[n-1]+h;
y[n_]:=y[n]=h f[x[n-1],y[n-1]]+y[n-1];
data=Table[{x[n],y[n]},{n,0,40}];
TableForm[data,TableHeadings->{None,{"(Subscript[x, n])","(Subscript[y, n])"}}]

Here we decrease the step size, and instead of outputting a table we output a plot.

Clear[x,y]
f[x_,y_]=(1+x)/(x*y);
x[0]=-3.5;
y[0]=0.1;
h = 0.001;
x[n_] := x[n] = x[n-1]+h;
y[n_]:= y[n] =  h f[x[n-1],y[n-1]]+y[n-1];
data = Table[{x[n],y[n]},{n,0,3800}];
eulerplot = ListPlot[data,BaseStyle->16,FrameLabel->{"x","y"},PlotStyle->Orange]

We can overlay the plot from Euler's Method against the direction field using the Show command.

Show[field,eulerplot]

Notice that the solution from Euler’s method seems to jump around near the end. This is because the solution does not exist where y=0 due to dy/dx=f(x,0) being infinite. So the derivative is infinite, and Euler’s method cannot capture that information. The points Euler’s method produces after this happens are meaningless.

DSolve to Solve DE

Mathematica has a built-in command to solve DEs. We have to tell Mathematica that y is a function of x when using DSolve. Note this is an equation we are solving, so there are double equal signs in the equation!

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

This produced two solutions, so let’s plot them. Notice as these are solutions to a DE they (both together) represent a family of curves. For the moment we will plot them separately to see why we got the two solutions instead of one. Also, notice Mathematica did not put in absolute values in the logarithms! I will make sure to do that, since if I don’t do that we will lose the part of the solution when x<0.

familya=-(Sqrt[2]) (Sqrt[x+C+Log[Abs[x]]])/.C->Range[-5,5,1/2]
familyb=(Sqrt[2]) (Sqrt[x+C+Log[Abs[x]]])/.C->Range[-5,5,1/2]

famplota=Plot[familya,{x,-4,4},BaseStyle->16,FrameLabel->{"x","y"},PlotStyle->{Red},Frame->True];
famplotb=Plot[familyb,{x,-4,4},BaseStyle->16,FrameLabel->{"x","y"},PlotStyle->{Blue},Frame->True];
Show[famplota,famplotb,PlotRange->All]

Both these blue and red curves together form the family of curves that solves the DE.

We can see whay they mean more clearly if we overlay them over the direction field.

Show[field,famplota,famplotb,PlotRange->All]

Notice that the arrows on the blue and red are “butting into” each other where the two curves meet.

Since the solution would follow the arrows, a solution cannot move from the red curve to the adjoining blue curve.

This is also precisely the place where the Euler's method became invalid.

Now that we understand why the solution was two sets of curves instead of one, we can combine them together to plot the single family of curves for the solution to the DE as an implicit function, (y^2)=2 (x+C+Log[Abs[x]] ). Plotting an implicit family of curves takes a bit more work, but here is one way you can do it.

family={};
Do[AppendTo[family,(y^2)==2 (x+C+Log[Abs[x]] )],{C,-5,5,1}]
family

famplot=ContourPlot[family,{x,-4,4},{y,-4,4},BaseStyle->16,FrameLabel->{"x","y"},ContourStyle->{Red},Frame->True]

Here is the family of curves overlaid with the direction field.

Show[field,famplot,PlotRange->All]

Here is the family of curves overlaid with the direction field and showing some streams.

Show[fieldstreams,famplot,PlotRange->All]

DSolve to Solve IVP

We can also use DSolve to solve an IVP

Clear[x,y]
DSolve[{y'[x]==(1+x)/(x*y[x]),y[-3.5]==0.1},y[x],x]

Mathematica is telling us with the error message that there is only one solution. Recall the solution to the DE had two solutions at this point.

This is directly related to what we saw earlier with the arrows in the direction field butting into each other at y=0.

There is a complex number in there ((i^2)=-1) but this is also inside a square root, so this function will output a real number for some x values! Notice for this to happen, I have to use Mathematica’s solution, with no absolute value inside the logarithm.

Here is a plot of the solution to the IVP, overlaid with the Euler’s solution we found earlier and the direction field. We can see the Euler method solution was pretty good.

IVPplot=Plot[(Sqrt[2]) (Sqrt[(2.252237031504632` -3.141592653589793` I)+x+Log[x]]),{x,-4,4},BaseStyle->16,FrameLabel->{"x","y"},PlotStyle->{Darker[Green]},Frame->True];

Show[field,eulerplot,IVPplot]

Show[field,eulerplot,IVPplot,PlotRange->{{-4,0},{0,4}}]