Lab 3 - Solving equations, tips and tricks

Some things we need to know when solving for stuff....

Andy Long, MAT229, Spring 2025

Exercise 0 -- Some important built-in functions

Special Constants

({{Standard notation, Mathematica notation}, {π≈3.14159, Pi}, {e≈2.71828, E}})

Important and Useful Commands

({{Functionality, Mathematica notation}, {plot the graph of a function, Plot[...]}, {square root of something, (Sqrt[…]), Sqrt[...]}, {absolute value of something, |…|, Abs[...]}, {sine of something (radian mode),sin(…), Sin[...]}, {cosine of something (radian mode),cos(…), Cos[...]}, {tangent of something (radian mode),tan(…), Tan[...]}, {natural logarithm of something, ln(…), Log[...]}, {inverse sine of something, (sin^(-1))(…), ArcSin[...]}, {inverse cosine of something, (cos^(-1))(…), ArcCos[...]}, {inverse tangent of something, (tan^(-1))(…), ArcTan[...]}, {integrate something, ∫(…)       
differentiate something, (d/(d□))(…)   , Integrate[...]
        D[...]}})

Exercise 1 -- Asymptotic Behavior

Featured:
	PolynomialQuotient
	assigning the output of an expression to a function with “=”, rather than “:=”

Your quiz yesterday featured a question about the asymptotic behavior of a rational function:
 h(x)=((3(x^3)-(x^2)+4x-3)/(x-1))
Let’s take a look at this function and how it approaches its asymptotic quadratic.... Plot the function and the asymptotic quadratic on the interval [-3,3].

The function “PolynomialQuotient” will accomplish the long division, and return the asymptotic function: I hope that it’s the same as the one you got on the quiz!

numerator[x_]:=3x^3- x^2+4x-3
denominator[x_]:=x-1
h[x_]:=numerator[x]/denominator[x]
asymptotic[x_]=PolynomialQuotient[numerator[x],denominator[x],x]
asymptotic[x]
Plot[{h[x],asymptotic[x]},{x,-4,4}]

Exercise 2 -- Solving for zeros....

Featured:
	Simplify
	Solve
	NSolve
	N -- converts an analytic expression to a decimal
	/. -- soliciting the value of the assignment to a variable in an expression
	% -- the previous expression

Oftentimes we need to find zeros of a derivative function, in order to find extrema of a function. Yesterday’s quiz featured a function like this:
 f(x)=4 x (e^(1-(x^2)))
(although I’ve switched the exponent on the exponential from x^2-1 to 1-x^2 for this exercise, because I need a function with some extrema!

Let’s find the derivative of f, call it fp, and plot the two functions on the interval [-3,3].

f[x_]:=4 x E^(-x^2+1)
f[x]
fp[x_]=f'[x] 
fp[x] (* Notice that the form of fp is cluttered -- not "simplified" So let's use Simplify to make it tidier, and to see clearly the zeros: *)
fp[x_]=Simplify[f'[x]]
fp[x]
plot=Plot[{f[x],fp[x]},{x,-3,3}]

Now we need to find the derivative of f, and the solve for where its derivative has zeros. In this case, it’s easy to do analytically. So we can check our answers with the basic solve commands for Mathematica, Solve and NSolve. 

-- Solve will solve the problem analytically, if possible, and provide an exact answer
-- NSolve will solve the problem numerically, and generally give approximate answers.

The general syntax requires double “==” for the equation, i.e. 
Solve [g[x]==0, x]
and note that we have to tell Mathematica which variable we’re solving for.

Solve[fp[x]==0,x]
NSolve[fp[x]==0,x]

Because there are two solutions, the answer is given as a list. Suppose that we want the first answer, and want to assign it to “min”. 

We will first create “solutions” to hold the list, and then ask for the first in the list.

solutions = Solve[fp[x]==0,x]
first = solutions[[1]]
min = x/.first
{min,f[min]}
N[%]

Exercise 3 -- Let’s find the inflection points, and add those and the extrema of the function to the plot

Featured:
	D -- to find higher derivatives

fpp[x_]=D[f[x],{x,2}]
solutions = Solve[fpp[x]==0,x]
second = solutions[[2]]
inflection = x/.second
{inflection,f[inflection]}
N[%]

Show[
plot,
ListPlot[
{{min,f[min]},{-min,f[-min]}},PlotStyle->{Black,PointSize[Large]}],
ListPlot[
{{0,0},{inflection,f[inflection]},{-inflection,f[-inflection]}},PlotStyle->{Red,PointSize[Large]}]
]

Exercise 4 -- Let’s see how calculus can help us to find the “best-fitting” line to some linear data.

Featured:
	ListPlot
	Transpose -- turns a pair of lists into a list of pairs
	Total -- adds up a list of data
	Mean -- computes the mean of a list of data
	Length -- the length of a list
	RandomReal, RandomVariate, NormalDistribution -- these are used to create 		random data of various types.

nRandomPoints=100;
xs=RandomReal[10,nRandomPoints]; (* 100 random real numbers between 0 and 10 *)
ys=Pi xs-E + RandomVariate[NormalDistribution[0,3],nRandomPoints];
 (* We're generating y values from the x values, and adding some random noise which is normally distributed *)
dataPlot=ListPlot[Transpose[{xs,ys}]]

Manipulate[
Show[
Plot[ a x+b,{x,0,10},PlotRange->{-10,40},PlotLabel->a x+b],dataPlot],
{a,0,5,.01},
{b,-10,10,.01}
]

Now we’re going to use calculus to find the optimal slope of the best fit line. We’re going to cheat a little, because there are really two “variables” in play here -- a and b -- and we’re doing univariate calculus -- just one variable. 

So I’m going to take the liberty of finding the slope of the best fit line that goes through the center of gravity of all of these points. For that reason, I’m going to define the center of gravity, as

n=Length[xs]
xmean = Total[xs]/nRandomPoints
ymean = Total[ys]/nRandomPoints
xmean = Mean[xs]
ymean = Mean[ys]

Now we define the “Objective function” that we’re going to minimize, to find the best fitting line. 

The linear model will be 

lm(x)=m(x-xmean) + ymean

The important unknown here is the slope m: we want to find that. We are going to use the data values, and compare the linear model’s predicted value, lm((Subscript[x, i])), to the actual data value (Subscript[y, i]). We’d like them to be the same, but they won’t be -- the data doesn’t actually follow a line. 

The function will be called 
objective(m) = (∑_(i=1)^n)(((lm((Subscript[x, i]))-(Subscript[y, i])))^2) 
	= (∑_(i=1)^n)(((m((Subscript[x, i])-xmean)+ymean-(Subscript[y, i])))^2)
This function will be smallest -- a minimum -- when m is chosen well. By computing the derivative of E, and setting it equal to zero, we’ll find our best fitting slope.

lm[x_]:=m(x-xmean) + ymean
objective[m_]:=Sum[(lm[xs[[i]]]-ys[[i]])^2,{i,1,nRandomPoints}]
Solve[D[objective[m],m]==0,m]
slope=m/.%[[1]]
bvalue=ymean-slope*xmean

Manipulate[
Show[
Plot[ a x+bvalue,{x,0,10},PlotRange->{-5,40},PlotLabel->slope x+bvalue],
ListPlot[{{xmean,ymean}},PlotStyle->PointSize[Large]],
dataPlot
],
{{a,slope},0,5,.01}
]