Algebraic Equations: Analytic Solutions using “Solve”

ClearAll["Global`*"]
(* First do some reminders of regular "Solve" *)
answers=Solve[p v == n R T,n];
n = n/.answers[[1]];
Print["n = ",n]

ClearAll["Global`*"]
answersa=Solve[x==x0 + v0x t,t]; (* solve horizontal projectile motion equation to find t in terms of x *)
t=t/.answersa[[1]];
answersb=Solve[{x==x0 + v0x t; y==y0 + v0y t- (1/2)g (t^2)},y]; (*use the prior answer in vertical equation to find y(x) instead of y(t) *)
y = y/.answersb[[1]]//Simplify;
Print["The trajectory equation is: y(x) = ",y]

Ordinary Differential Equations: Analytic Solutions using “DSolve”

	Suppose I want to work on a certain differential equation, such as y’ + 3xy = 2. It turns out that this problem has an analytical solution, which means that we really shouldn’t be using Mathematica for this problem to start with. However, suppose that, as a result of alien influences, you’ve lost your memory of how to do kindergarten math. Let’s see if we can find some ways to get Mathematica to help us. It turns out that this is Mathematica’s strong suit.
	The above equation is a first order equation (only y’ appears), so we’ll need one initial or boundary condition before we can get started. 
	For my problem, the initial condition will be given as y(0) = 1. 
	An answer will consist of a description for y(x), either as a list of numbers for x with the corresponding values of y, or, in the best case, an actual function for y(x).
	First, let’s look at the DSolve function to try to get an analytical solution. The first argument of DSolve is a {list} containing the original equation along with all of the boundary conditions. Each of these equations must use “==” in place of a mere “=”. We can just use a regular “prime” to indicate the derivative. Note that the prime character (an apostrophe) appears between the name of the function and it’s dependence. So, you type y’[x],             but NOT y[x]’                 and NOT y’.
	The second argument of DSolve is the name of the output function, with dependencies listed explicitly, and the final argument is the name of the independent variable. 
	The output of the DSolve function is a “rule”, so we’ll have to use the slash-dot notation to convert it into a “real” Mathematica equation using a second line. For the example, recall that we want y’ + 3xy = 2, with y(0) = 1.

ClearAll["Global`*"]
yrule=DSolve[{y'[x]+3 x y[x] == 2, y[0]==1},y[x],x];
y1[x_]=y[x]/.yrule[[1]];
Print["y(x) = ",y1[x]]
Plot[y1[x],{x,0,8},PlotRange-> {{0,8},{0,1.5}}]

In this plot, I only bothered to plot between x = 0 and x = 8, but since it's a symbolic solution, it will actually be valid for *any* value of x (like, x = 9, for example).
	You’ll notice that the analytic solution uses a new function called “Erfi”. By definition, Erfi[x] = Erf[x i]/i. Since we expect that Erf[some imaginary number] is itself imaginary, then Erf[x i] divided by i should end up being a real number. So, when x is real, then we expect Erfi[x] will also be real. For example:

Print["Erfi(3i) = ",Erfi[3.0 I]];Print["Erf(3i) = ",Erf[3.0 I]]
Print["Erf(3) = ",Erf[3.0]];Print["Erfi(3) = ",Erfi[3.0]]
Print["Erfi(0.2i) = ",Erfi[0.2 I]];Print["Erf(0.2) = ",Erf[0.2]]

Ordinary Differential Equations: Numeric Solutions using “NDSolve”

OK, so we now have an “analytical” (or algebraic) solution. Analytic solutions are ALWAYS more true than numeric solutions, when they exist. But sometimes, problems don’t even have an analytical solution. For such problems, we  want to use the function NDSolve, where the "N" obviously stands for “Numerical”. For such problems, we have to specify ahead of time the range of the independent variable (x, in this case) that we want to obtain solutions for. When Mathematica gives us an answer, it will only be valid for the range specified. So, I’ll now use NDSolve to solve the above problem for x again, this time for x between 0 and 8, as plotted above.  Whatever answer comes out will NOT be useable at x = 9! 
	When you use NDSolve, the answers that Mathematica  generates are actually just two lists of numbers: the first is for x (in this case the list ranges from 0 to 8), and a list of numbers for y that correspond to the values of x in the first list. Mathematica automatically does two nice things when making these lists:
	First, it chooses values of x that are close enough to each other to “smoothly” represent the answer, as needed, based on the shape of the answer. So, for places where the result has a higher slope on the plot, Mathematica has more data points for the answer.
	Second, the way it generates the lists is such that when you ask for y[x] for a value of x that isn’t actually on the list, Mathematica will automatically look at the nearest available values of x and interpolate between them without even bothering you with the details of how it’s doing it. As we’ll see soon, this can have some bad side effects if you are not paying attention to what you’re doing!

yrule2=NDSolve[{y2'[x]+3 x y2[x] -2 ==0, y2[0]==1},y2[x],{x,0,8}];

y2[x_]=y2[x]/.yrule2[[1]];

I should probably check some of the answers to make sure that y1 came out the same as y2:

Print["The symbolic value of y(0.5) was ",y1[0.5]," and the numerical value was ",y2[0.5]]
Print["The symbolic value of y(2) was ",y1[2.0]," and the numerical value was ",y2[2.0]]

Well, that looks encouraging. We can also plot the results to see them visually:

plot1=Plot[y1[x],{x,0,8},PlotStyle->Directive[RGBColor[1,0,0],AbsoluteThickness[6]]];
plot2=Plot[y2[x],{x,0,8},PlotStyle->Directive[RGBColor[0,0,1],AbsoluteThickness[2]]];
Legended[Show[plot1,plot2,AxesStyle->Black,FrameLabel->{"x","y"},Frame->True,FrameTicks->All,PlotRange-> {{0,8},{0,1.5}},BaseStyle->{12,FontFamily->"Helvetica"}],Placed[LineLegend[{Red,Blue,Darker[Green],Orange},{Style["Symbolic",12],Style["Numeric",12]},LegendMargins->{{15,15},{5,5}},LabelStyle->Directive[{Bold}],LegendFunction->Framed],{{0.575,0.5},{0,0}}]]

What would happen if we tried to find y2[20] (or, some other x value outside the range we specified earlier)?

Print["The symbolic value of y(9) was ",y1[9.0]," and the numerical value was ",y2[9.0]]
Print["The symbolic value of y(20) was ",y1[20.0]," and the numerical value was ",y2[20.0]]

That’s interesting! It knows it can’t do this problem, but it tried to guess an answer based on some of the answers that it did have. As you can see, Mathematica did a really bad job guessing... the correct answer for y(when x=20) was 0.03336, but it guessed 0.151. Thank goodness we have the error message that tells us not to trust this answer! Let’s try to see this visually...

plot3=Plot[y1[x],{x,0,20},PlotStyle->Directive[RGBColor[1,0,0],AbsoluteThickness[6]],PlotRange-> {{0,40},{0,1.5}}];
plot4=Plot[y2[x],{x,0,20},PlotStyle->Directive[RGBColor[0,0,1],AbsoluteThickness[2]],PlotRange-> {{0,40},{0,1.5}}];
Legended[Show[plot3,plot4,ImageSize->{600,300},AxesStyle->Black,FrameLabel->{"x","y"},Frame->True,FrameTicks->All,PlotRange-> {{0,20},{0,1.5}},BaseStyle->{12,FontFamily->"Helvetica"}],Placed[LineLegend[{Red,Blue,Darker[Green],Orange},{Style["Symbolic = True solution",12],Style["Numeric... extended range",12]},LegendMargins->{{15,15},{5,5}},LabelStyle->Directive[{Bold}],LegendFunction->Framed],{{0.375,0.5},{0,0}}]]

From looking at this plot, it appears that even though the numeric solution can’t really be trusted for x>8, it looks like the numeric answers aren’t really that bad until x of about 12 or 13. But they just get worse and worse after that!

Second Order Differential Equations

Let’s do a word problem. A mass m = 5 kg oscillates horizontally on a spring having k = 10 N/m, and experiences some friction that is proportional to speed having a constant of proportionality b = 1 Ns/m. The initial conditions are x(0) = 1 m, and v(0) = 0 m/s. In other words, we start by pulling the mass 1 m to the side and then let go. 
It isn’t even clear yet that we need a differential equation. Let’s write Newton’s Second Law for the mass:
ΣF = ma
-kx - bv = ma
ma + bv + kx = 0
mx’’ + b x’ + kx = 0
OK, now that’s clearly a differential equation. Since the highest derivative is second order, we need two initial (or boundary) conditions, and that’s how many we have, so we’re good. Let’s try to get both an analytical solution (I’ll call it x1), and a numerical one (x2). As a side note, this is EXACTLY the kind of problem you should already be able to do quickly WITHOUT Mathematica!

(* get the symbolic solution with DSolve ... *)
xrule1=DSolve[{m x1''[t]+b x1'[t] +k x1[t] ==0, x1[0]==1,x1'[0]==0},x1[t],t];
x1[t_]=x1[t]/.xrule1[[1]]//FullSimplify

Of course, we’ll need *numbers* for every parameter to eventually perform an NDSolve. Let’s use these:

m=5; k = 10; b = 1.0;x0=1;v0=0;
Print["With numbers, x1(t) = ",x1[t]//FullSimplify]

I note that somehow, these numbers were such that the hyperbolic functions got converted into “regular” functions!  Now let’s do the “numeric” version of the answer.

xrule2=NDSolve[{m x2''[t]+b x2'[t] +k x2[t] ==0, x2[0]==1,x2'[0]==0},x2[t],{t,0,20}];
x2[t_]=x2[t]/.xrule2[[1]]; (* this is the numeric solution with NDSolve ... *)

plot1=Plot[x1[t],{t,0,20},PlotStyle->Directive[RGBColor[1,0,0],AbsoluteThickness[6]]];
plot2=Plot[x2[t],{t,0,20},PlotStyle->Directive[RGBColor[0,0,1],AbsoluteThickness[2]]];
Legended[Show[plot1,plot2,ImageSize->{600,300},AxesStyle->Black,FrameLabel->{"Time (s)","Position (m)"},Frame->True,FrameTicks->All,PlotRange-> {{0,20},{-1.25,1.25}},BaseStyle->{12,FontFamily->"Helvetica"}],Placed[LineLegend[{Red,Blue,Darker[Green],Orange},{Style["Symbolic",12],Style["Numeric",12]},LegendMargins->{{15,15},{5,5}},LabelStyle->Directive[{Bold}],LegendFunction->Framed],{{0.575,0.025},{0,0}}]]

We could again examine what would happen if you tried to look at points that are outside the range specified in our NDSolve command:

plot1=Plot[x1[t],{t,0,30},PlotStyle->Directive[RGBColor[1,0,0],AbsoluteThickness[6]]];
plot2=Plot[x2[t],{t,0,30},PlotStyle->Directive[RGBColor[0,0,1],AbsoluteThickness[2]]];
Legended[Show[plot1,plot2,ImageSize->{600,300},AxesStyle->Black,FrameLabel->{"Time (s)","Position (m)"},Frame->True,FrameTicks->All,PlotRange-> {{0,30},{-1.25,1.25}},BaseStyle->{12,FontFamily->"Helvetica"}],Placed[LineLegend[{Red,Blue,Darker[Green],Orange},{Style["Symbolic = True",12],Style["Numeric",12]},LegendMargins->{{15,15},{5,5}},LabelStyle->Directive[{Bold}],LegendFunction->Framed],{{0.55,0.05},{0,0}}]]

Example: A Relativistic ODE that we’ve already solved (Prob 8-1-7)

For this problem, a mass is accelerated by a constant force. It starts at rest at the origin. Let’s use m = 0.5 kg and F = 30 N; both are pretty modest values. We’ll let it compute results for about a year. The “ClearAll” is there so that Mathematica doesn’t become confused about what x, y, or t mean, since we’ve used them all already above.

ClearAll["Global`*"]
m = 0.5;F = 30.0; c = 2.998 (10^8);tmax = 3 (10^7);

You’ll recall that last week, most of the problem was spent discovering the equation itself, as opposed to solving it. Here’s what we found for the equation; I’ve also just copied our Analytical solution (and I called it x1) so we have something to compare against. Now, I’ll also do this problem numerically:

x1[t_]=((c/F)) (Sqrt[(m^2)(c^2)+(F^2)(t^2) ]) - ((m (c^2))/F); (* copied previous solution we found by hand *)
xrule=NDSolve[{((m/(((1-(((((x'[t])/c)))^2)))^(3/2))))x''[t]-F ==0, x[0]==0,x'[0]==0},x[t],{t,0,tmax}];
x2[t_]=x[t]/.xrule[[1]];

plot1=Plot[x1[t],{t,0,tmax},PlotStyle->Directive[RGBColor[1,0,0],AbsoluteThickness[6]]];
plot2=Plot[x2[t],{t,0,tmax},PlotStyle->Directive[RGBColor[0,0,1],AbsoluteThickness[2]]];
Legended[Show[plot1,plot2,ImageSize->{600,300},AxesStyle->Black,FrameLabel->{"Time (s)","Position (m)"},Frame->True,FrameTicks->All,BaseStyle->{12,FontFamily->"Helvetica"}],Placed[LineLegend[{Red,Blue,Darker[Green],Orange},{Style["Symbolic",12],Style["Numeric",12]},LegendMargins->{{15,15},{5,5}},LabelStyle->Directive[{Bold}],LegendFunction->Framed],{{0.55,0.05},{0,0}}]]

It’s clear that  Mathematica’s numerical answer agrees with our “paper” or “hand-generated” answer pretty well, but as with all numeric solutions, it’s too bad it doesn’t give a nice symbolic answer. As a check, we expect that the slope of this plot is c when evaluated at large times (the mass accelerates continuously, slowly approaching the speed of light):

Print["Speed at large time for hand solution is ",x1'[tmax]]
Print["Speed at large time for numeric solution is ",x2'[tmax]]

This is great! Who would have guessed that Mathematica could give the derivative of a list of numbers as easily as it does for a true function? Clearly, we haven’t quite reached the speed of light yet after about one year, but we’re already pretty close. 

Now, let’s evaluate the quality of the numerical solution in a more quantitative way.

error[t_] = ((x2[t]-x1[t])/(x1[t])); (* this is relative error. "error" of 0.1 means 10% error, and so on... we can only do this because we have a true solution, x1  *)
Plot[error[t],{t,0,tmax},PlotStyle->Directive[RGBColor[1,0,0],AbsoluteThickness[2]],ImageSize->{800,300},AspectRatio->Full]

That’s pretty impressive... the results are accurate to about 5 parts in one billion (99.9999995% accuracy). 
We used our “hand generated solution” for this comparison, but of course we could use DSolve (instead of NDSolve) to get that solution, too. 
To do this, we’ll have to reset all the numerical values using another ClearAll:

ClearAll["Global`*"]
xrule=DSolve[{((m/(((1-(((((x'[t])/c)))^2)))^(3/2))))x''[t]-F ==0, x[0]==0,x'[0]==0},x[t],t]
x[t_]=Simplify[x[t]/.xrule[[2]]]//Expand

You should carefully note the [[2]] in the above line... it turns out that this problem has multiple analytical solutions: the first one is negative, and the second is positive. We want the positive solution!

You might recall that when we solved this problem by hand, we did so by first finding v[t]. Here, you see both the answer we got for v[t], as well as the original ODE we solved to find it:

ClearAll["Global`*"]
m = 0.5;F = 30.0; c = 2.998 (10^8);tmax = 3 (10^7);
v1[t_]= (Sqrt[((((F c t))^2)/((((F t))^2)+(((m c))^2))) ]);
vrule=NDSolve[{((m/(((1-(((((v[t])/c)))^2)))^(3/2))))v'[t]-F ==0, v[0]==0},v[t],{t,0,tmax}];
v2[t_]=v[t]/.vrule[[1]];

plot1=Plot[v1[t],{t,0,tmax},PlotStyle->Directive[RGBColor[1,0,0],AbsoluteThickness[6]],PlotRange-> {0,c}];
plot2=Plot[v2[t],{t,0,tmax},PlotStyle->Directive[RGBColor[0,0,1],AbsoluteThickness[2]],PlotRange-> {0,c}];
Legended[Show[plot1,plot2,ImageSize->{600,300},AxesStyle->Black,FrameLabel->{"Time (s)","Velocity (m/s)"},Frame->True,FrameTicks->All,BaseStyle->{12,FontFamily->"Helvetica"}],Placed[LineLegend[{Red,Blue,Darker[Green],Orange},{Style["Symbolic",12],Style["Numeric",12]},LegendMargins->{{15,15},{5,5}},LabelStyle->Directive[{Bold}],LegendFunction->Framed],{{0.55,0.05},{0,0}}]]

It’s easy to see again that the speed is asymptotically approaching the speed of light.

Example: An ODE that we’ve already solved (Prob 8-1-7)
This next problem is the “spherical” heat problem we did a few days ago in class. 
You can’t use DSolve directly, because this a first order equation with 2 boundary conditions; the second BC is needed because some of the physical constants weren’t actually given!

ClearAll["Global`*"]
Trule=DSolve[{Q2== 4 π (r^2) T'[r],T[Ra]== Ta},T[r],r]; (* Q2 is like (Qdot/k) *)
T[r_]=T[r]/.Trule[[1]];Print["T(r) = ",T[r]]
(* we now have T(r) to within a constant "Q2". Next, let's find Q2...*)
qrule=Solve[T[Rb]==Tb,Q2];Q2=Q2/.qrule[[1]];
Print["Q = ",Q2];Print["T(r) without Q2 in it is ",T[r]//Simplify]
Ra=1;Rb=2;Ta=100.0;Tb=0;
Print["With numbers put in, T(r) is ",T[r]//Simplify]

ClearAll["Global`*"] (* this is just to see what goes wrong if you try to do this kind of strange problem  all in one step... *)
Ra=1;Rb=2;Ta=100.0;Tb=0;
Trule=DSolve[{Q2== 4 π (r^2) T'[r],T[Ra]== Ta, T[Rb]==Tb},T[r],r]

ClearAll["Global`*"] (* But of course, some problems can be done all in one step... *)
y[x_]=y[x]/.DSolve[{y''[x]+3y'[x]==6,y[0]==2,y'[0]==0},y[x],x][[1]]