The mathematics 4 course in Wolfram Mathematica

Preface

All students of our Faculty have to take the courses of mathematics. To go through this process, they use books, university textbooks, materials from lectures, etc. We wanted to help them by creating some interactive study material. For this reason, the project processes one of the courses of mathematics of our Faculty ‒ the Mathematics 4 course.
The project is divided into four main parts. The first part deals with the particular topics of Mathematics 4 course. For each topic, the most important definitions and theorems are given and some specific examples solved. The source code written in Mathematica is included.
In the next part, interactive visualizations of some notions and problems are developed. They could be used not only by students but also by their lecturers as an interactive tool during classes. They are developed in the spirit of demonstrations at [] but tailored to the students of our Faculty.
As there is no built-in function for calculating bivariate limits in Mathematica, the third part of the project contains a function which can calculate some of them.
Finally, as a practical application of the skills developed during the Mathematics 4 course, a physical problem is solved ‒ namely the problem of ballistic curves.
As has been already said, the main goal of the project is to create an interactive tool for students of our Faculty. They can see that there is another way how to solve a problem and that Mathematica could help them with the task they are not able to solve due to difficult calculations. Furthermore, some more difficult and interesting problems are solved.
The whole project will be available on the Internet and is written as a notebook in Mathematica in order that students could work with it and could change the content in any way they want.

Introduction to Wolfram Mathematica

. Wolfram Mathematica

Wolfram Mathematica, or simply Mathematica, is a computational software program based on symbolic mathematics. It was invented by Stephen Wolfram in 1988 and is developed by Wolfram Research. The latest version is available for three major platforms ‒ Microsoft Windows, MacOS X, and Linux (supporting both 32-bit and 64-bit implementations) [].
Mathematica is divided into two main parts. The code is interpreted by the kernel, whereas the front end provides the graphical user interface. Much of the content can be edited interactively or generated by an algorithm [].
Programming in Mathematica differs from programming in other languages. The code consists of primitives. For a new user of Mathematica, it can be a problem to find their right combination which would result in the simplest solution. Furthermore, Mathematica has changed a lot from the first version even though the core principles remain the same. Obsolete functions have been replaced or discarded. As many highly specialized functions have been added, the built-in library has increased, and even the most difficult tasks can be very often resolved in a very simple way [].
Mathematica is not used only in mathematics and physics. Due to its strength, user-friendly interface, and possibility of creating interactive modules, it is often utilized in other fields such as biology or chemistry. Moreover, the source code can be accompanied by the text and the document can be formatted with ease, so it is possible to write whole articles using Mathematica only. Later versions also enable connection with other applications through MathLink protocol, parallel processing, dynamic interactivity, and high-performance computing. This makes Mathematica a complex tool with a wide range of applications [].

. WolframAlpha

Another project developed by Wolfram Research is a computational knowledge engine WolframAlpha. Unlike search engines, it does not give a list of websites containing the searched information but tries to give the answer directly by searching the web and using some dynamic computations [].
WolframAlpha is based on Mathematica (its language and built-in functions). Additional data is gathered from various websites. For instance, some reports about users of Facebook are generated. They contain how often they use some specific word, how popular their friends are, etc. [].
Since WolframAlpha is available on the Internet free of charge, it is very popular among students.

The Mathematics 4 course in examples

. Differential equations

The Mathematics 4 course deals with three types of ordinary differential equations:

In contrast to the hand calculation, Mathematica can handle the above mentioned types of differential equations with ease. Just use the command DSolve to find the solution.

Example ..: Find the solution of the equation y'-2xy=3(x^2)(e^(x^2))passing through the point M=(0,5) [].

The command DSolve is used to search the solution. The initial condition is represented by the second equation.

sol=DSolve[{y'[x]-2x*y[x]==3(x^2)*(E^(x^2)),y[0]==5},y[x],x]

Let us have a look what the function y looks like.

Plot[y[x]/.sol[[1]],{x,-2,1},PlotStyle->Black,ImageSize->300]

The domain of y is the set of all real numbers.

Clear[sol];

Example ..: Find the solution of the equation y'=cos(x) passing through the point M=(0,1).

Again, the command DSolve is used.

DSolve[{y'[x]==Cos[x],y[0]==1},y[x],x]

The domain of y are all real numbers.

Example ..: Find the general solution of the equation y''+2y'+y=3(e^(-x))(Sqrt[x+1]) [].

DSolve[y''[x]+2*y'[x]+y[x]==3*(E^(-x))*(Sqrt[x+1]),y[x],x]

As it can be seen from the solution, Dom(y)=(-1,+∞).

. Multivariable limits

Definition ..:
Let M⊂(ℝ^n) and (OverVector[a])∈(ℝ^n). We say that (OverVector[a]) is a limit point of M if and only if (Subscript[U, ϵ]^( * ))((OverVector[a]))⋂M≠∅ for all ϵ >0, where (Subscript[U, ϵ]^( * ))((OverVector[a])) is the punctured ϵ-neighbourhood of (OverVector[a]).

Definition ..:
Let f:(ℝ^n)->ℝ and (OverVector[a]) be a limit point of the domain of f. The limit of f, as (OverVector[x]) tends to (OverVector[a]), is  c∈ℝ, written as (lim_((Overscript[x, ->])->(Overscript[ a, ->])))f((OverVector[x]))=c if and only if for every ϵ>0 there exists δ>0 such that for all (OverVector[x]) ∈(Subscript[D, f]),
0<‖(OverVector[x])-(OverVector[a])‖<δ implies |f((OverVector[x]))-c|<ϵ.

Definition ..:
Let f:(ℝ^n)->ℝ, M⊂(ℝ^n), and (Overscript[a, ->]) be a limit point of (Subscript[D, f])⋂M. The limit of f, as (OverVector[x]) tends to (OverVector[a]), with respect to the set M is c∈ℝ, written as (lim_((OverVector[x])-> (OverVector[a]), (OverVector[x])∈M )) f((OverVector[x]))=c, if and only if for every ϵ>0 there exists δ>0 such that for all (OverVector[x])∈(Subscript[U, δ]^( * ))((OverVector[a]))⋂(Subscript[D, f])⋂M we have|f((OverVector[x]))-c|<ϵ.

Theorem ..:
Let f:(ℝ^2)->ℝ, (OverVector[a])=((Subscript[a, x]),(Subscript[a, y])). If the limit of f, as (OverVector[x]) tends to (OverVector[a]), is equal to c∈ℝ and the sequential limits exist and are finite, then

(lim_(x -> (Subscript[a, x])))((lim_(y -> (Subscript[a, y])))f((OverVector[x])))=(lim_(y -> (Subscript[a, y])))((lim_(x -> (Subscript[a, x])))f((OverVector[x])))=c.

Corollary ..:
If both sequential limits exist and are finite but different, the overall limit does not exist.

Theorem ..:
Let f:(ℝ^n)->ℝ, (OverVector[a]) be a limit point of (Subscript[D, f]), and (lim_((OverVector[x])->(OverVector[a])))f((OverVector[x])) exists. Then for any A⊂(ℝ^n) such that (OverVector[a]) is a limit point of A,

(lim_((OverVector[x])->(OverVector[a]), (OverVector[x])∈A))f((OverVector[x]))=(lim_((OverVector[x])->(OverVector[a])))f((OverVector[x])) .

Corollary ..:
Suppose A⊂(ℝ^n), B⊂(ℝ^n), and (OverVector[a]) is a limit point of them. If

(lim_((OverVector[x])->(OverVector[a]), (OverVector[x])∈A))f((OverVector[x]))≠(lim_((OverVector[x])->(OverVector[a]), (OverVector[x])∈B))f((OverVector[x])) ,

Example ..: Let g(x,y)=((xy^2)/((x^2)+(y^4))). Show that (lim_((x,y)->(0,0)))g(x,y) does not exists.

g[x_,y_]=((x*(y^2))/((x^2)+(y^4)));

The non-existence of the limit could be confirmed using sequential limits. If the sequential limits existed and were different, the limit of g would not exist.

Limit[Limit[g[x,y],y->0],x->0]

Limit[Limit[g[x,y],x->0],y->0]

The sequential limits are the same. Therefore, the decision whether the limit exists cannot be made using sequential limits. Nonetheless, the limit would be equal to zero provided it existed. Another method could be used to prove the non-existence of a limit. For instance, approaching (0, 0) via points which lie on the lines x=ky or y=kx.

Limit[g[k*y,y],y->0]

Limit[g[x,k*x],x->0]

Again, the result provides no information on the non-existence of the limit. Finally, when we choose parabolae x=(ky^2), we will get a significant result.

Limit[g[k*(y^2),y],y->0]

The numerical value of the limits with respect to different parabolae clearly depends on k. Consequently, the overall limit does not exist. The level curves of the function g have the shape of parabolae x=(ky^2). Along these parabolae, f attains different constant values.

Clear[g];

Example ..: Let f(x,y)=((2xy)/(xy+2x-y)). Show that (lim_((x,y)->(0,0)))f(x,y) does not exists [].

f[x_,y_]=((2*x*y)/(x*y+2*x-y));

At first, let us have a look what the sequential limits look like.

Limit[Limit[f[x,y],y->0],x->0]

Limit[Limit[f[x,y],x->0],y->0]

Since they are the same, another method have to be used to prove the non-existence. Try to approach the point (0, 0) via lines y=kx.

Limit[f[x,y]/.y->k*x,x->0]

Although the result is that the limit is 0 on these lines, it is not always true.

f[x,y]/.y->k*x//Simplify

It is possible to use L'Hospital's rule for k=2.

((D[2 k*x,x])/(D[2+k*(-1+x),x]))

To summarize, the limit does not exist because we were able to find two different paths on which the limit is not the same.

Clear[f];

. Partial derivative

Definition ..:
Suppose f:(ℝ^n)->ℝ is defined in a neighbourhood of (OverVector[a])=((Subscript[a, 1]),(Subscript[a, 2]),...,(Subscript[a, n])). Partial derivative of f at the point (OverVector[a]) with respect to the i-th variable is defined as

((∂f)/(∂(Subscript[x, i])))((OverVector[a])):=(lim_(h->0))((f((Subscript[a, 1]),(Subscript[a, 2]),...,(Subscript[a, i-1]),(Subscript[a, i])+h,(Subscript[a, i+1]),...,(Subscript[a, n]))-f((Subscript[a, 1]),(Subscript[a, 2]),...,(Subscript[a, i-1]),(Subscript[a, i]),(Subscript[a, i+1]),...,(Subscript[a, n])))/h)

Definition ..:
By a gradient of a function f at the point (OverVector[a]) is understood

grad f((OverVector[a])):=(((∂f)/(∂(Subscript[x, 1])))((OverVector[a])),...,((∂f)/(∂(Subscript[x, n])))((OverVector[a])))

Remark ..:
The symbol ∇ f is used for the gradient of f.

Note: The partial derivative of f at the point (OverVector[a]) with respect to x has a geometric interpretation of the slope of the tangent at the point ((OverVector[a]),f((OverVector[a]))) to the curve given by the intersection of the plot of f with the plane that is perpendicular to the domain of f and contains x-axis. Similarly, the partial derivative of f at (OverVector[a]) with respect to y. These facts can be nicely illustrated on the function f=(Sqrt[1-(x^2)-(y^2)]) at some point.

f[x_,y_]=(Sqrt[1-(x^2)-(y^2)]);

a={0.5,0.4};

a3D=Append[a,f@@a];

gradient=Grad[f[x,y],{x,y}]/.{x->a3D[[1]],y->a3D[[2]]}

hemisphere=Plot3D[f[x,y],{x,-1,3},{y,-1,3},
Mesh->{{a3D[[1]]},{a3D[[2]]}},MeshStyle->Thickness[0.003],
FaceGrids->{{0,0,-1}},FaceGridsStyle->Directive[Dotted],
BoxRatios->{4,4,1},AxesLabel->{x,y,z},
ColorFunction->"Aquamarine"];

The tangents are defined because we know the starting point and the direction vector.

line1[x_]=gradient[[1]]*(x-a3D[[1]])+a3D[[3]];

line2[y_]=gradient[[2]]*(y-a3D[[2]])+a3D[[3]];

tangents={ParametricPlot3D[{x,a3D[[2]],line1[x]},{x,-1,3},PlotStyle->{Thickness[0.003],Black}],
ParametricPlot3D[{a3D[[1]],y,line2[y]},{y,-1,3},PlotStyle->{Thickness[0.003],Black}]};

Other objects are defined to make the illustration more lucid.

planes=Graphics3D[{FaceForm[GrayLevel[0.5]],EdgeForm[],Opacity[0.3],
Polygon[{{-1,a3D[[2]],0},{-1,a3D[[2]],1},{3,a3D[[2]],1},{3,a3D[[2]],0}}],
Polygon[{{a3D[[1]],3,0},{a3D[[1]],3,1},{a3D[[1]],-1,1},{a3D[[1]],-1,0}}]}];

lines=Graphics3D[{Black,Dashed,Thickness[0.003],
Line[{{-1,a3D[[2]],0},{3,a3D[[2]],0}}],
Line[{{a3D[[1]],-1,0},{a3D[[1]],3,0}}]}];

point1=x/.Solve[line1[x]==0];

point2=y/.Solve[line2[y] ==0];

angle1=ParametricPlot3D[{{x,a3D[[2]],(Sqrt[0.4-(((x-point1))^2)])}},{x,1,1.15},PlotStyle->{Thickness[0.003],Black}];

angle2=ParametricPlot3D[{{a3D[[1]],y,(Sqrt[0.4-(((y-point2))^2)])}},{y,1.2,1.31},PlotStyle->{Thickness[0.003],Black}];

text1=Graphics3D[{Style[Text["(Subscript[α, x])",{1.2,a3D[[2]],0.15}],18],Style[Text["(Subscript[t, x])",{1,a3D[[2]],0.7}],18]}];

text2=Graphics3D[{Style[Text["(Subscript[α, y])",{a3D[[1]],1.5,0.13}],18],Style[Text["(Subscript[t, y])",{a3D[[1]],1.2,0.6}],18]}];

Show[{hemisphere, planes,tangents,lines,angle1,angle2,text1,text2},ViewPoint->{4,3,1.3},ImageSize->500]

Show[{hemisphere,planes, tangents,lines,angle1,text1},ViewPoint->{0,∞,0},Ticks->None,AxesLabel->{x,None,z}]

Show[{hemisphere,planes, tangents,lines,angle2,text2},ViewPoint->{∞,0,0},Ticks->None,AxesLabel->{None,y,z}]

Clear[f,a2D,a,gradient,hemisphere,line1,line2,tangents,planes,lines,point1,point2,angle1,angle2,text1,text2];

Example ..: Calculate ((∂f)/(∂x)) and ((∂f)/(∂y)) of the function f=(x^3)(y^((y^4)+y)).

f[x_,y_]=(x^3)*(y^((y^4)+y));

parcDerX=Simplify[D[f[x,y],x]]

parcDerY=Simplify[D[f[x,y],y]]

Clear[f,parcDerX,parcDerY];

. Directional derivative

Definition ..:
Suppose f:(ℝ^n)->ℝ is defined in a neighbourhood of (OverVector[a]) and (OverVector[s]) is a non-zero vector in (ℝ^n). A directional derivative of f along (OverVector[s]) at (OverVector[a]) is defined as

((∂f)/(∂(OverVector[s])))((OverVector[a])):=(1/(‖(OverVector[s])‖))(lim_(h->0))((f((OverVector[a]) + h(OverVector[s]))-f((OverVector[a])))/h)

Theorem ..:
Suppose f:(ℝ^n)->ℝ and all partial derivatives at (OverVector[a]) are continuous. Then f has derivatives in all directions and

((∂f)/(∂(OverVector[s])))((OverVector[a]))=(1/(‖(OverVector[s])‖))⟨grad f((OverVector[a])),(OverVector[s])⟩.

Example ..: Find the directional derivative of f=arctan(xy) along (OverVector[s])=((Sqrt[2]),(Sqrt[2])) at the point (OverVector[a])=(0,1).

f[x_,y_]=ArcTan[x*y];

s={(Sqrt[2]),(Sqrt[2])};

a={0,1};

Definition of the function which can calculate directional derivative of the function of two variables follows:

directionalDerivative[f_,s_,a_]:=(Grad[f[x,y],{x,y}] /.{x->a[[1]],y->a[[2]]}).s/Norm[s];

directionalDerivative[f,s,a]

Clear[f,s,a];

Example ..: Find the directional derivative of f=({{{, ({{0, x=-y}, {(x+y)sin((1/((x+y)))), x≠-y}})}})  along (OverVector[s]) =(1,-1)  at the point (OverVector[a])=(0,0).

f[x_,y_]=If[x==-y,0,(x+y)*Sin[(1/((x+y)))]];

s={1,-1};

a={0,0};

Plot3D[f[x,y],{x,-0.1,0.1},{y,-0.1,0.1},ColorFunction->"Aquamarine",ViewPoint-> {0,1,1},ImageSize->300]

Equation () can be used if partial derivatives of f at the point (OverVector[a]) with respect to x and y exist.

derivativeX=D[f[x,y],x]/.{x->a[[1]],y->a[[2]]}

derivativeY=D[f[x,y],y]/.{x->a[[1]],y->a[[2]]}

Mathematica gave us the result that ((∂f)/(∂x))((OverVector[a]))=0 and ((∂f)/(∂y))((OverVector[a]))=0. Nonetheless, it can be proved that it is incorrect by using the definition of a partial derivative.

derX=Limit[(1/h)((f[x,y]/.{x->a[[1]]+h,y->a[[2]]})-(f[x,y]/.{x->a[[1]],y->a[[2]]})),h->0]

derY=Limit[(1/h)((f[x,y]/.{x->a[[1]],y->a[[2]]+h})-(f[x,y]/.{x->a[[1]],y->a[[2]]})),h->0]

Mathematica returns an interval because the limit does not exist. Thus, the directional derivative must be calculated directly from the definition.

dirDer=(1/(Norm[s]))*Limit[
(1/h)((f[x,y]/.{x->a[[1]]+h*s[[1]],y->a[[2]]+h*s[[2]]})-(f[x,y]/.{x->a[[1]],y->a[[2]]})),h->0]

Clear[f,s,a,derivativeX,derivativeY,derX,derY,dirDer];

Example ..: Find the gradient of g=((x^3)+x+y)(e^(-(x^4)-(y^2)))-2(-(x^5)-(1/5)(y^2)+x)(e^(-(x^2)-(y^2))). Plot g with several gradient vectors. Find directional derivatives at (OverVector[a])=(-1.4,0.4) along several directions.

g[x_,y_]=((x^3)+x+y)*(E^(-(x^4)-(y^2)))-2*(-(x^5)-(1/5)*(y^2)+x)*(E^(-(x^2)-(y^2)));

a={-1.4,0.4};

gradient=Simplify[Grad[g[x,y],{x,y}]];

gradientAtA=gradient /.{x->a[[1]],y->a[[2]]};

Plot3D[{g[x,y]},{x,-2.5,2.5},{y,-2.5,2.5}, ColorFunction->"Aquamarine",AxesLabel->{x,y,z},MeshStyle->Darker[Blue],ImageSize->250,ViewPoint->{-2,-2.5,1}]

contours1=ContourPlot[g[x,y],{x,-2.5,2.5},{y,-2.5,2.5},Contours->14, ColorFunction->"Aquamarine"];

contours2=ContourPlot[g[x,y],{x,-2.5,0.5},{y,-1.5,1.5},Contours->11, ColorFunction->"Aquamarine"];

gradientField=VectorPlot[gradient,{x,-2.5,2.5},{y,-2.5,2.5},VectorPoints->15,VectorStyle->Arrowheads[0.02]];

Show[{contours1,gradientField},ImageSize->250]

The gradient of g is perpendicular to level curves at all points. Its length depends on the growth of the function at the specific point.

scaleParameter=1.5;

arrow=Graphics[Line[{{-0.5,1/4},{0,0},{-0.5,-1/4}}]];

As we demonstrated in the previous figure, the gradient is perpendicular to the level curves of the function. We can also show that the gradient at (OverVector[a]) is greater than all directional derivatives at the same point.

gradientEndPoint={gradientAtA[[1]]/scaleParameter+a[[1]],gradientAtA[[1]]/scaleParameter+a[[2]]};

gradientGraphics=Graphics[{Red,Arrowheads[{{.05,1,arrow}}],Arrow[{a,gradientEndPoint}]}];

angles=Table[2*π*(n/12),{n,12}];

directions=Table[{Cos[angles[[n]]],Sin[angles[[n]]]},{n,12}];

dirDerivatives=Table[directionalDerivative[g,directions[[n]],a],{n,12}];

dirDerPoints1=Table[{Cos[angles[[n]]]*dirDerivatives[[n]]/scaleParameter+a[[1]],Sin[angles[[n]]]*dirDerivatives[[n]]/scaleParameter+a[[2]]},{n,6}];

dirDerPoints2=Table[{(-1)*Cos[angles[[n]]]*dirDerivatives[[n]]/scaleParameter+a[[1]],(-1)*Sin[angles[[n]]]*dirDerivatives[[n]]/scaleParameter+a[[2]]},{n,6}];

dirDerGraphics={Graphics[{{Arrowheads[{{.05,1,arrow}}],
Table[Arrow[{a,dirDerPoints1[[n]]}],{n,6}]}}],
Graphics[{{Arrowheads[{{.05,0,arrow}}],
Table[Arrow[{dirDerPoints2[[n]],a}],{n,6}]}}]};

Labeled[Show[{contours2,gradientGraphics,dirDerGraphics},ImageSize->250],PointLegend[{Red,Black},{"gradient","directional derivatives"},LegendMarkers->"⟶" ],Right]

The graph shows directional derivatives of g at (OverVector[a]) along several directions. The length of each arrow corresponds to the numerical value of the directional derivative along the arrow. The greatest growth of the function is in the direction of the gradient.

Clear[g,a,gradient,gradientAtA,contours1,contours2,gradientField,scaleParameter,arrow,gradientEndPoint,gradientGraphics,angles,directions,dirDerivatives,dirDerPoints1,dirDerPoints2,dirDerGraphics];

. Local extrema

Definition ..:
Suppose f:(ℝ^n)->ℝ is defined on M⊂ (ℝ^n) and (OverVector[a])∈M. Then f attains a local minimum or a local maximum at (OverVector[a]) with respect to M if and only if there is U((OverVector[a])) such that for all (OverVector[x])∈M⋂U((OverVector[a])), f((OverVector[x]))≥f( (OverVector[a]))or f((Overscript[x, ->]))≤f( (OverVector[a])), respectively.
Similarly, f attains a sharp local minimum or a sharp local maximum at (OverVector[a]) if and only if for all (OverVector[x])∈M⋂(SuperStar[U])((OverVector[a])), f((OverVector[x]))>f( (OverVector[a])) or f((OverVector[x]))<f( (OverVector[a])), respectively.
A point which is either a local minimum or a local maximum is called a local extremum.

Theorem .. (The first derivative test):
Let (OverVector[a]) be a local extremum of f and all partial derivatives of f at (OverVector[a]) exist. Then

((∂f)/(∂(Subscript[x, i])))((OverVector[a]))=0

Definition ..:
Suppose f:(ℝ^n)->ℝ and all the second partial derivatives of f at (OverVector[a]) exist and are continuous. A Hessian matrix of f at (OverVector[a]) is defined as a square matrix of second partial derivatives

H(Subscript[(f), ij])((OverVector[a])):= (((∂^2)f)/(∂(Subscript[x, i])∂(Subscript[x, j])))((OverVector[a])).

Theorem .. (The second derivative test):
Let H(f)((OverVector[a])) be the Hessian matrix of f at (OverVector[a]) and (OverVector[a]) be a critical point.
1. If H(f)((OverVector[a])) is a positive definite matrix, then (OverVector[a]) is a local minimum point.
2. If H(f)((OverVector[a])) is a negative definite matrix, then (OverVector[a]) is a local maximum point.
3. If H(f)((OverVector[a])) is indefinite matrix, then (OverVector[a]) is a saddle point.

Remark ..:
1. H(f)((OverVector[a])) is a positive definite matrix if and only if all its eigenvalues are positive.
2. H(f)((OverVector[a])) is a negative definite matrix if and only if all its eigenvalues are negative.
3. H(f)((OverVector[a])) is indefinite matrix if and only if some of its eigenvalues are positive and some negative.

Example ..: Find the local extrema and the saddle points of f=((x^2)-(y^2))(e^((-(x^2)-(y^2))/2)) with respect to M={(x,y)∈(ℝ^2):(x^2)+(y^2)<9}.

f[x_,y_]=((x^2)-(y^2))*(E^((-(x^2)-(y^2))/2));

The local extremum points can be roughly estimated from the figure.

Plot3D[{f[x,y]},{x,-3,3},{y,-3,3}, ColorFunction->"BlueGreenYellow",ImageSize->250,AxesLabel->{x,y,z},ViewPoint->{1,-1.5,1.5}]

The first derivative test is used in order to find the critical points.

a=Solve[D[f[x,y],x]==0 &&D[f[x,y],y]==0 && (x^2)+(y^2)<9,{x,y}]

We found five critical points ((Subscript[(OverVector[a]), 1]), ..., (Subscript[(OverVector[a]), 5])). Lets us try to use the second derivative test. We calculate Hessian matrices and analyze whether (Subscript[(OverVector[a]), i]), i ∈ {1, ..., 5} is a local extremum point or not.

hessianMatrix[f_,a_]:={{D[f[x,y],x,x]/.a,D[f[x,y],x,y]/.a},{D[f[x,y],x,y]/.a,D[f[x,y],y,y]/.a}}

secondDerivativeTest[mat_]:=Module[{eigenvs=Eigenvalues[mat]},If[And@@Positive[eigenvs],Print["Local minimum."],If[And@@Negative[eigenvs],Print["Local maximum."],If[And@@NonPositive[eigenvs],Print["The second derivative test cannot be used."],If[And@@NonNegative[eigenvs],Print["The second derivative test cannot be used."],Print["Saddle point."]]]]]];

secondDerivativeTest[hessianMatrix[f,a[[1]]]]

secondDerivativeTest[hessianMatrix[f,a[[2]]]]

secondDerivativeTest[hessianMatrix[f,a[[3]]]]

secondDerivativeTest[hessianMatrix[f,a[[4]]]]

secondDerivativeTest[hessianMatrix[f,a[[5]]]]

Now the result can be demonstrated graphically using the ContourPlot command.

contours=ContourPlot[f[x,y],{x,-3,3},{y,-3,3},Contours->10,ColorFunction->"BlueGreenYellow",PlotLegends->Automatic,AxesLabel->{x,y}];

criticalPoints=Graphics[Table[{Text[Style[(Subscript["a", i]),Black,17],{(x/.a)[[i]],(y/.a)[[i]]+0.25}],PointSize[Large],Black,Point[{(x/.a)[[i]],(y/.a)[[i]]}]},{i,5}]];

Labeled[Show[{contours ,criticalPoints},ImageSize->250],PointLegend[{Black},{"critical points"},LegendMarkers->{Graphics[Point[{0,0}]]}],Right]

There are also built-in functions FindMinimum and FindMaximum in Mathematica which search for local extrema numerically. However, these functions return only one local extremum point and sometimes are not able to find any solution.

Clear[f,a,contours,criticalPoints,hessianMatrix];

Example ..: Find the local extrema and the saddle points of f=(x^3)+(y^2)+(1/2)(z^2)-3xz-2y+2z [].

f[x_,y_,z_]=(x^3)+(y^2)+(1/2)*(z^2)-3*x*z-2*y+2*z;

We will use the same procedure as in the example ...

a=Solve[D[f[x,y,z],x]==0 &&D[f[x,y,z],y]==0 &&D[f[x,y,z],z]==0 ,{x,y,z}]

hessianMatrix[f_,a_]:={{D[f[x,y,z],x,x]/.a,D[f[x,y,z],x,y]/.a,D[f[x,y,z],x,z]/.a},{D[f[x,y,z],x,y]/.a,D[f[x,y,z],y,y]/.a,D[f[x,y,z],y,z]/.a},
{D[f[x,y,z],z,x]/.a,D[f[x,y,z],y,z]/.a,D[f[x,y,z],z,z]/.a}}

The function secondDerivativeTest from the previous example can be used again.

secondDerivativeTest[hessianMatrix[f,a[[1]]]]

secondDerivativeTest[hessianMatrix[f,a[[2]]]]

Clear[f,a,hessianMatrix,secondDerivativeTest];

. Constrained extrema

Definition ..:
Suppose f:(ℝ^n)->ℝ is defined on M⊂ (ℝ^n) and (OverVector[a])∈M. Then f attains a local minimum or a local maximum at (OverVector[a]) with respect to M if and only if there is U((OverVector[a])) such that for all (OverVector[x])∈M⋂U((OverVector[a])), f((OverVector[x]))≥f( (OverVector[a]))or f((Overscript[x, ->]))≤f( (OverVector[a])), respectively.
Similarly, f attains a sharp local minimum or a sharp local maximum at (OverVector[a]) if and only if for all (OverVector[x])∈M⋂(SuperStar[U])((OverVector[a])), f((OverVector[x]))>f( (OverVector[a])) or f((OverVector[x]))<f( (OverVector[a])), respectively.
A point which is either a local minimum or a local maximum is called a local extremum.

Theorem .. (The first derivative test):
Suppose f,(Subscript[g, 1]),(Subscript[g, 2]),..., (Subscript[g, m]) are functions of n variables and have continuous partial derivatives. Suppose that f has a local extremum with constraints (Subscript[g, s])=0, s∈{1,...,m}, at the point (OverVector[a]) and (Subscript[▽g, s])((OverVector[a]))≠0 for all s∈{1,...,m}. Then there are numbers (Subscript[λ, 1]),(Subscript[λ, 2]),...,(Subscript[λ, m]) such that

▽f((OverVector[a]))=(Subscript[λ, 1]) (Subscript[▽g, 1])((OverVector[a]))+(Subscript[λ, 2]) (Subscript[▽g, 2])((OverVector[a]))+...+(Subscript[λ, m]) (Subscript[▽g, m])((OverVector[a])).

Remark ..:
It is convenient to rewrite () as follows,

▽L((OverVector[a]))=0,

Example ..: Find the local extrema of  f=(1/x)+(1/y) subjected to the constraint (1/(x^2))+(1/(y^2))-(1/4)=0 [].

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

g[x_,y_]=(1/(x^2))+(1/(y^2))-(1/4);

At first, we try the built-in functions Minimize and Maximize.

Minimize[{f[x,y],g[x,y]==0},{x,y}]

Maximize[{f[x,y],g[x,y]==0},{x,y}]

It can be seen that these functions work just fine in our example. Nevertheless, let us try to solve the same problem mimicking a hand calculation. The first derivative test helps us to find the critical points together with corresponding values of the Lagrange multiplier.

a=Solve[D[f[x,y]-λ*g[x,y],x]==0 &&D[f[x,y]-λ*g[x,y],y]==0 &&g[x,y]==0 ,{x,y,λ}]

In the next step, the Hessian matrix of the Lagrange function is constructed at the critical points.

hessianMatrix[f_,a_]:={{D[f,x,x]/.a,D[f,x,y]/.a},{D[f,x,y]/.a,D[f,y,y]/.a}};

Finally, the function secondDerivativeTest is defined to identify the type of an extremum from the eigenvalues of the Hessian matrix.

secondDerivativeTest[mat_]:=Module[{eigenvs=Eigenvalues[mat]},
If[And@@Positive[eigenvs],Print["Local minimum."],
If[And@@Negative[eigenvs],Print["Local maximum."],
Print["Our basic second derivative test is not sufficient. One has to reduce the Hessian matrix to the subspace that is tangent to the constraints and then use the second derivative test on the resulting matrix that is effectively lower-dimensional."]]]];

secondDerivativeTest[hessianMatrix[f[x,y]-λ*g[x,y],a[[1]]]]

secondDerivativeTest[hessianMatrix[f[x,y]-λ*g[x,y],a[[2]]]]

Next, it will be graphically demonstrated that at the points where local constrained extrema are attained, the gradient of f is just a scalar multiple of the gradient of g, i.e., () holds. To this purpose, denote the critical points as follows:

a1={x/.a[[1]],y/.a[[1]]};

a2={x/.a[[2]],y/.a[[2]]};

Moreover, introduce a function that returns the gradient vector of a function f at a point startPoint shifted to this startPoint.

arrows[f_,startPoint_,scaleParameter_,colour_]:=
Module[{gradient,endPoint},
gradient=Grad[f[x,y],{x,y}]/.{x->startPoint[[1]],y->startPoint[[2]]};
endPoint={gradient[[1]]*scaleParameter+startPoint[[1]],gradient[[2]]*scaleParameter+startPoint[[2]]};
Graphics[{colour,Arrow[{startPoint,endPoint}]}]]

Finally, define all needed graphical objects and combine all of them into an illustrative figure.

contoursf=ContourPlot[f[x,y],{x,-5,5},{y,-5,5},Contours->{-2,-1,-(1/(Sqrt[2])),-0.5,0,0.5,(1/(Sqrt[2])),1,2},ContourShading->False,ContourStyle->Blue,MaxRecursion->3];

contoursg=ContourPlot[g[x,y]==0,{x,-5,5},{y,-5,5},ContourShading->False,ContourStyle->Darker[Green]];

b={-2 (Sqrt[2]),2 (Sqrt[2])};

text=Graphics[{Table[{Text[Style[(Subscript["a", i]),Black,12],{(x/.a)[[i]]+0.1,(y/.a)[[i]]-0.5}],PointSize[Large],Point[{(x/.a)[[i]],(y/.a)[[i]]}]},{i,2}],
Text[Style["b",Black,12],{b[[1]],b[[2]]-0.5}],PointSize[Large],Point[{b[[1]],b[[2]]}],
Text[Style["▽f((Subscript[a, 1]))",Blue,10],{a1[[1]]-1.2,a1[[2]]-0.5}],
Text[Style["▽g((Subscript[a, 1]))",Darker[Green],10],{a1[[1]],a1[[2]]+0.7}],
Text[Style["▽f((Subscript[a, 2]))",Blue,10],{a2[[1]]-0.6,a2[[2]]-1.5}],
Text[Style["▽g((Subscript[a, 2]))",Darker[Green],10],{a2[[1]]-0.65,a2[[2]]-0.02}],
Text[Style["▽f(b)",Blue,10],{b[[1]]-0.3,b[[2]]-1}],
Text[Style["▽g(b)",Darker[Green],10],{b[[1]]+1,b[[2]]-0.1}],
}];

Show[{contoursf ,contoursg,text,arrows[f,a1,10,Blue],arrows[g,a1,10,Darker[Green]],arrows[f,a2,10,Blue],arrows[g,a2,10,Darker[Green]],arrows[f,b,10,Blue],arrows[g,b,10,Darker[Green]]},ImageSize->250]

Gradients of f and g at (Subscript[a, 1]) (local minimum point) are proportional. The same holds true at (Subscript[a, 2]) (local maximum point). On the other hand, this is not true at the point b which is not a critical point.

Clear[f,g,a,a1,a2,b,contoursf,contoursg,text];

Example ..: Find the local extrema of f(x,y)=((x^2)-(y^2))(e^((-(x^2)-(y^2))/2)) subjected to the constraint (((x-1))^2)+(((y-1))^2)=(1/2).

f[x_,y_]=((x^2)-(y^2))*(E^((-(x^2)-(y^2))/2));

g[x_,y_]= (((x-1))^2)+(((y-1))^2)-(1/2);

It can be seen that the constraint is the equation of the circle with the following parameters:

ax=1;

ay=1;

radius=(Sqrt[(1/2)]);

Again, let us start with built-in function.

Minimize[{f[x,y],g[x,y]==0},{x,y}]

The function f is too complicated and so is its gradient.Therefore, it is not possible to solve the equation for the critical points symbolically. So try the numerical version of the above used functions.

NMinimize[{f[x,y],g[x,y]==0},{x,y}]

NMaximize[{f[x,y],g[x,y]==0},{x,y}]

This works pretty well, both for the minimum and the maximum. Nevertheless, we could also try to use the first derivative test.

a=NSolve[g[x,y]==0&&D[f[x,y]-λ*g[x,y],x]==0 &&D[f[x,y]-λ*g[x,y],y]==0 ,{x,y,λ},Reals]

Again, we use the functions hessianMatrix and secondDerivativeTest from the previous example in order to find the Hessian matrix of the Lagrange function at the critical points and identify the type of an extremum.

secondDerivativeTest[hessianMatrix[f[x,y]-λ*g[x,y],a[[1]]]]

secondDerivativeTest[hessianMatrix[f[x,y]-λ*g[x,y],a[[1]]]]

Clear[f,g,ax,ay,radius,a,arrows,hessianMatrix,secondDerivativeTest];

Example ..: The Baraboo, Wisconsin, plant of International Widget Co. uses aluminum, iron, and magnesium to produce high-quality widgets. The quantity of widgets which may be produced using x tons of aluminum, y tons of iron, and z tons of magnesium is Q(x,y,z)=xyz. The cost of raw materials is: aluminum, $6 per ton; iron, $4 per ton; and magnesium, $8 per ton. How many tons of aluminum, iron, and magnesium should be used to manufacture 1000 widgets at the lowest possible cost? []

Considering the cost of raw materials, the function f can be defined.

f[x_,y_,z_]=6*x+4*y+8*z;

q[x_,y_,z_]=x*y*z-1000;

We want to minimize the cost, i.e., the function f. Q is the constraint. The Lagrange function can be defined as:

lagrange[x_,y_,z_]=f[x,y,z]-λ*q[x,y,z];

The first derivative test is used to find the solution. The fourth equation is given by the constraint. The result is approximately:

a=Solve[D[lagrange[x,y,z],x]==0 && D[lagrange[x,y,z],y]==0 
&& D[lagrange[x,y,z],z]==0 &&q[x,y,z]==0
,{x,y,z,λ},Reals]//N

Clear[f,q,lagrange,a];

. Double integral

In this section, we will consider only double integrals. The transition to higher dimensions can be made in a similar way.

Definition ..:
A finite sequence a=(Subscript[x, 0])<(Subscript[x, 1])<...<(Subscript[x, n])=b, where a<b, is called a partition of the interval ⟨a,b⟩.

Definition ..:
Let {(Subscript[x, 0]),(Subscript[x, 1]),...,(Subscript[x, n])} be a partition of the interval ⟨a,b⟩ and {(Subscript[y, 0]),(Subscript[y, 1]),...,(Subscript[y, m])} be a partition of the interval ⟨c,d⟩. Then (Subscript[{(Subscript[S, i,j]):=⟨(Subscript[x, i-1]),(Subscript[x, i])⟩ ×⟨(Subscript[y, i-1]),(Subscript[y, i])⟩ }, i=1,j=1]^(n,m)) is called a partition of the interval ⟨a,b⟩×⟨c,d⟩. Additionally, measure of the interval I=⟨a,b⟩×⟨c,d⟩ is defined as μ(I)=(b-a)(d-c).

Definition ..:
Suppose f:(ℝ^2)->ℝ is bounded on the interval I=⟨a,b⟩×⟨c,d⟩ and D is a partition of this interval. Let (Subscript[v, i,j])=(inf_(x∈(Subscript[S, i,j]))) f(x) and (Subscript[V, i,j])=(sup_(x∈(Subscript[S, i,j])))f(x). Then L(D,f):=(∑_(i=1,j=1)^(n,m))(Subscript[v, i,j])μ((Subscript[S, i,j])) is called a lower sum for f on D. Similarly, U(D,f):=(∑_(i=1,j=1)^(n,m))(Subscript[V, i,j])μ((Subscript[S, i,j])) is called an upper sum for f on D.

Definition ..:
Suppose f:(ℝ^2)->ℝ is bounded on the interval I=⟨a,b⟩×⟨c,d⟩, M⊂(ℝ^2), and D is a partition of this interval. Then (sup_M)L(D,f) or (inf_M) U(D,f) is called a lower Riemann integral or an upper Riemann integral, respectively. If (sup_M) L(D,f) =(inf_M) U(D,f), the function f is said to be Riemann integrable over I and (∫_I)f=(sup_M) L(D,f).

Definition ..:
Let I=⟨a,b⟩×⟨c,d⟩ and D ={(Subscript[S, ij])} be a partition of this interval. Then a norm of D is defined as ‖D‖:=max μ((Subscript[S, ij])), where  0≤i≤n and  0≤j≤m .

Theorem ..:
Suppose f:(ℝ^2)->ℝ is continuous on the interval I=⟨a,b⟩×⟨c,d⟩ Then (∫_I)f exists and if (Subscript[{(Subscript[D, n])}, n=1]^∞) is a partition sequence such that (lim_(n->∞))‖(Subscript[D, n])‖=0, then

(∫_I)f =(lim_(n->∞))L((Subscript[D, n]),f)=(lim_(n->∞))U((Subscript[D, n]),f).

Theorem .. (The Fubini theorem):
Suppose f:(ℝ^2)->ℝ is continuous on the interval I=⟨a,b⟩×⟨c,d⟩ Then

(∫_I)f(x,y)dxdy=(∫_a^( b))((∫_c^( d))f(x,y)dy)dx=(∫_c^( d))((∫_a^( b))f(x,y)dx)dy.

Definition ..:
Let (OverVector[F]):(ℝ^n)->(ℝ^s) and (OverVector[a]) be a point of its domain such that for all i∈{1,...,s} and j∈{1,...,n} there exists ((∂(Subscript[F, i]))/(∂(Subscript[x, j])))((OverVector[a])). Then Jacobian matrix is, by definition,

((∂(OverVector[F]))/(∂(OverVector[x])))((OverVector[a]))≡((D((Subscript[F, 1]),...,(Subscript[F, s])))/(D((Subscript[x, 1]),...,(Subscript[x, n])))):=(({{((∂(Subscript[F, 1]))/(∂(Subscript[x, 1]))), ⋯, ((∂(Subscript[F, 1]))/(∂(Subscript[x, n])))}, {⋮, ⋱ ,  }, {((∂(Subscript[F, s]))/(∂(Subscript[x, 1]))),  , ((∂(Subscript[F, s]))/(∂(Subscript[x, n])))}}))((OverVector[a])).

Definition ..:
A function (OverVector[F]):(ℝ^n)->(ℝ^n) is called regular over an open set M∈(ℝ^n) if and only if it has a continuous Jacobian matrix and det ((D((Subscript[ϕ, 1]),...,(Subscript[ϕ, s])))/(D((Subscript[ξ, 1]),...,(Subscript[ξ, n]))))((OverVector[x]))≠0 for all (OverVector[x])∈M.

Theorem .. (The substitution theorem):
Suppose (OverVector[ϕ]):(ℝ^n)->(ℝ^n), (OverVector[ϕ]):((Subscript[ξ, 1]),...,(Subscript[ξ, n]))->((Subscript[x, 1]),...,(Subscript[x, n])) is an injective and regular vector function from P⊂(ℝ^n) to Q⊂(ℝ^n). Then for arbitrary set M⊂Q

(∫_M)f((Subscript[x, 1]),...,(Subscript[x, n]))d(Subscript[x, 1]),...,d(Subscript[x, n])=(∫_((ϕ^(-1))(M)))f((OverVector[ϕ])((Subscript[ξ, 1]),...,(Subscript[ξ, n]))) |det ((D((Subscript[ϕ, 1]),...,(Subscript[ϕ, s])))/(D((Subscript[ξ, 1]),...,(Subscript[ξ, n]))))|d(Subscript[ξ, 1]),...,d(Subscript[ξ, n]).

Note: In the substitution theorem, the absolute value of Jacobian determinant appears. Substitutions to the polar, spherical, and cylindrical coordinates systems are used very often. Therefore, here are Jacobian determinants of them:

JacobianDeterminant[f_List,vars_List]:=
If[Equal@@Map[Dimensions,{f,vars}],Simplify[Det[Outer[D,f,vars]]],"The dimensions of input parameters are not the same."]

polar=JacobianDeterminant[{ρ*Cos[ϕ],ρ*Sin[ϕ]},{ρ,ϕ}]

spherical=JacobianDeterminant[{ρ*Cos[θ]*Cos[ϕ],ρ*Cos[θ]*Sin[ϕ],ρ*Sin[θ]},{ρ,θ,ϕ}]

As θ runs from -(π/2)to (π/2), cos(θ) is non-negative. Thus, the absolute value of the determinant is (ρ^2) cos(θ) for every ρ and θ.

cylindrical=JacobianDeterminant[{ρ*Cos[ϕ],ρ*Sin[ϕ],h},{ρ,ϕ,h}]

Clear[ϕ,θ,h,polar,spherical,cylindrical];

Example ..: Find the surface area bounded by the curves y=2(x^2) and y=1-x.

f[x_]=2*(x^2);
g[x_]=1-x;

To find the limits of integration, we need to locate intersections of the plots of f and g.

sol=Solve[f[x]==g[x],x]

Show[{Graphics[{Style[Text[x/.sol[[1]],{(x/.sol[[1]])+0.25,-0.25}],12],Style[Text[x/.sol[[2]],{(x/.sol[[2]])+0.25,-0.25}],12]}],Plot[{f[x]/.a->1,g[x]/.a->1},{x,-1.1,1.1},PlotStyle->{Blue,Green},PlotLegends->{"f[x]","g[x]"}]},Axes->True,Ticks->None,GridLinesStyle->Dashed,ImageSize->200,GridLines->{{x/.sol[[1]],x/.sol[[2]]},{}}]

The limits of integration are -1 and (1/2) (with respect to x) and the functions f and g (with respect to y). Since the surface area is searched, the functions which is to be integrated is just 1.

Integrate[Integrate[1,{y, 2*(x^2),1-x}],{x,-1,(1/2)}]

This example could be also solved without any knowledge of double integrals. We can calculate the integral of both g and f from -1 to (1/2), and the result would be the difference between them.

Clear[f,g,sol];

Example ..: Calculate the integral (∫_( ))(∫_M) (1/(Sqrt[2x-(x^2)]))dxdy, where M is an area bounded by the curves x=0, (y^2)=-x+2.

f[x_]=(1/(Sqrt[2*x-(x^2)]));

Let us have a look at M to find the limits of integration:

Show[{Graphics[{Style[Text["2",{2.5,-0.5}],12],Style[Text["0",{0.5,-0.5}],12]}],
ContourPlot[{x==0,(y^2)==-x+2},{x, -3,3},{y,-3,3},ContourStyle->{Blue,Green},PlotLegends->{"x = 0","(y^2) = - x + 2"}]},Axes->{True,False},Frame->False,Ticks->None,ImageSize->200,GridLinesStyle->Dashed,GridLines->{{2},{}}]

Considering the x-axis symmetry, only the area above the x-axis can be integrated and then multiplied by 2.

2*Integrate[Integrate[f[x],{y,0,(Sqrt[-x+2])}],{x,0,2}]

Clear[f];

Example ..: Find the integral ∫(∫_M) ln((x^2)+(y^2))dxdy, where M is an area in the first quadrant bounded by the curves (x^2)+(y^2)=1 and (x^2)+(y^2)=4 [].

f[x_,y_]=Log [(x^2)+(y^2)];

Labeled[Show[{Graphics[{Style[Text["0",{0.25,-0.25}],12],
Style[Text["1",{1.25,-0.25}],12],Style[Text["2",{2.25,-0.25}],12]}],
Graphics[{Green,Disk[{0,0},2,{0,Pi/2}],White,Disk[{0,0},1,{0,Pi/2}],Blue,Circle[{0,0}],Circle[{0,0},2]}]},
Axes->True,Frame->False,Ticks->None,ImageSize->200],PointLegend[{Green},{"integrated area"},LegendMarkers->"■"],Right]

Due to the polar symmetry, it is convenient to use the polar coordinates via the substitution x=ρ cos(ϕ) and x=ρ sin(ϕ).

fPolar[ρ_,ϕ_]=TransformedField["Cartesian"->"Polar",f[x,y],{x,y}->{ρ,ϕ}]//Simplify

The integration is much easier now. As it can be deduced from the figure, ρ runs from 1 to 2 and ϕ from 0 to 2π. The Jacobian determinant for polar substitution is ρ.

Integrate[Integrate[fPolar[ρ,ϕ]*JacobianDeterminant[{ρ*Cos[ϕ],ρ*Sin[ϕ]},{ρ,ϕ}],{ρ,1,2}],{ϕ,0,(π/2)}]//Simplify

Clear[f,fPolar];

Example ..: Find the integral (∫_( ))(∫_( ))(∫_B) (((1+(x^2)+(y^2)+(z^2)))^(-(1/2)))dxdydz, where B is an sphere with radius 1 and centre at the origin [].

g[x_,y_,z_]=(((1+(x^2)+(y^2)+(z^2)))^(-(1/2)));

The integration domain is a sphere so we use spherical coordinates system.

gSpherical[ρ_,θ_,ϕ_]=TransformedField["Cartesian"->"Spherical",g[x,y,z],{x,y,z}->{ρ,θ,ϕ}]//Simplify

The absolute value of Jacobian determinant is, in this case, (ρ^2)cos(θ). Since the integration domain is a sphere of radius 1, ρ goes from 0 to 1, ϕ from 0 to 2π, and θ from -(π/2)to (π/2) .

Integrate[Integrate[Integrate[gSpherical[ρ,θ,ϕ]*
Abs[JacobianDeterminant[{ρ*Cos[θ]*Cos[ϕ],ρ*Cos[θ]*Sin[ϕ],ρ*Sin[θ]},{ρ,θ,ϕ}]],{ρ,0,1}],{ϕ,0,2*π}],{θ,-(π/2),(π/2)}]

Clear[g,gSpherical];

Example ..: Find the integral (∫_( ))(∫_( ))(∫_M) (ze^((x^2)+(y^2)))dxdydz, where M={(x,y,z)∈(ℝ^3):(x^2)+(y^2)≤4,2≤z≤3} [].

f[x_,y_,z_]=z*(E^((x^2)+(y^2)));

We convert the function f to the cylindrical coordinates.

fCylindrical[ρ_,ϕ_,h_]=TransformedField["Cartesian"->"Cylindrical",f[x,y,z],{x,y,z}->{ρ,ϕ,h}]//Simplify

The integration domain M is a cylinder with radius 2. Thus, ρ runs from 0 to 2. Its height h goes from 2 to 3 and ϕ must run around full circle. The Jacobian determinant is ρ.

Integrate[Integrate[Integrate[fCylindrical[ρ,ϕ,h]*
JacobianDeterminant[{ρ*Cos[ϕ],ρ*Sin[ϕ],h},{ρ,ϕ,h}],{ρ,0,2}],{ϕ,0,2*π}],{h,2,3}]//Simplify

Clear[f,fCylindrical,JacobianDeterminant];

Interactive visualizations

. Solution to the differential equation

The demonstration depicts the plot of the solution of a chosen differential equation passing through a point (OverVector[a]) that is determined by dragging the locator. The equations are of the type  y'=f(x,y), where f:(ℝ^2)->ℝ is a continuous function. The vector field in the background represents the direction of the tangent line to the plot of the solution at each point. This field is given just by the right-hand side of the equation. The particular solution is fixed by the initial condition.

Manipulate[
equation=Switch[function,
"y",y'[x]==y[x],
"cos(x)",y'[x]==Cos[x],
"3(x^2)(E^(x^2))+2xy",y'[x]==3*(x^2)*(E^(x^2))+2x*y[x]
];
fun=Switch[function,
"y",y,
"cos(x)",Cos[x],
"3(x^2)(E^(x^2))+2xy",3*(x^2)*(E^(x^2))+2x*y
];
sol=DSolve[{equation,y[a[[1]]]==a[[2]]},y[x],x];
Show[{Plot[y[x]/.sol[[1]],{x,-5,5},PlotRange->5,PlotStyle->Black],
VectorPlot[(({1,fun})/(Sqrt[1+(fun^2)])),{x,-5,5},{y,-5,5},Axes->True,Frame->None,VectorPoints->15,VectorScale->Small,VectorStyle->{Opacity[0.3],Arrowheads[0.02],Thickness[0.001]}]},ImageSize->220,PlotRange->5,AspectRatio->1],
Row[{Control[{{a,{1,1}},Locator}],
Row[{Labeled[Control[{{function,"y",""},{"y","cos(x)","3(x^2)(E^(x^2))+2xy"}}],{"y' ="},{Left}],Spacer[60],Pane[Row[{" (OverVector[a]) = (",Dynamic[Round[a[[1]],0.1]],", ",Dynamic[Round[a[[2]],0.1]],")"}],ImageSize->80]}]}],ControlType->PopupMenu]

Clear[equation,fun,sol,a,function,y];

. Directional derivative

This applet demonstrates geometric interpretation of the directional derivative of the function f along a given vector (OverVector[s]) at the point (OverVector[a]). The upper inset in the figure shows the graph of f; the point (OverVector[a]) is depicted in blue, and the blue line lies in the direction of (OverVector[s]). The intersection of the graph with the plane, which is perpendicular to the xy-plane and contains the blue line, is displayed as a dashed line. The lower inset shows a two-dimensional cut together with the tangent to this intersection. The value of the directional derivative of f along (OverVector[s]) at (OverVector[a]) is just the tangent function of the angle between the tangent and the blue line. As one varies (OverVector[s]) (by choosing the value of ϕ, which denotes the angle between (OverVector[s]) and the x-axis), the directional derivative changes correspondingly.

Manipulate[Module[{f,a,a3D,directionalDerivative,plotf,s,horizontalLine,tangent,intersections,intersection1,intersection2,cut,plane,point,tangent2,cut2,plane2,point2},
f[x_,y_]=(Sqrt[(1-(1/2)*(x^2)-(y^2))]);
a={0.5,0.7};
a3D=Append[a,f@@a];
s={Cos[ϕ],Sin[ϕ]};
directionalDerivative[f_,s_,a_]:=(Grad[f[x,y],{x,y}] /.{x->a[[1]],y->a[[2]]}).s;
plotf=Plot3D[f[x,y],{x,-2,3},{y,-1,2},
PlotRange->{0,1},Mesh->None,AxesLabel->Automatic,
FaceGrids->{{0,0,-1}},FaceGridsStyle->Directive[Dotted],BoxRatios->{4,4,1},
ColorFunction->"Aquamarine",PlotStyle->Opacity[0.7]];
horizontalLine=ParametricPlot3D[{a3D[[1]]+s[[1]]*t,a3D[[2]]+s[[2]]*t,0},
{t,-10,10},PlotStyle->{Thickness[0.003],Blue}];
tangent=ParametricPlot3D[{a3D[[1]]+s[[1]]*t,a3D[[2]]+s[[2]]*t,a3D[[3]]+directionalDerivative[f,s,a3D]*t},
{t,-10,10},PlotStyle->{Thickness[0.003],Black}];
intersections=t/.NSolve[f[a3D[[1]]+s[[1]]*t,a3D[[2]]+s[[2]]*t]==0,t];
intersection1={a3D[[1]]+s[[1]]*intersections[[1]],a3D[[2]]+s[[2]]*intersections[[1]]};
intersection2={a3D[[1]]+s[[1]]*intersections[[2]],a3D[[2]]+s[[2]]*intersections[[2]]};
cut=ParametricPlot3D[{a3D[[1]]+s[[1]]*t, a3D[[2]]+s[[2]]*t,f[a3D[[1]]+s[[1]]*t,a3D[[2]]+s[[2]]*t]},
{t,intersections[[1]],intersections[[2]]},PlotStyle->{Thickness[0.003],Dashed}];
plane=Graphics3D[{FaceForm[GrayLevel[0.5]],EdgeForm[],Opacity[0.5],
Polygon[{{intersection1[[1]],intersection1[[2]],0},{intersection1[[1]],intersection1[[2]],1},
{intersection2[[1]],intersection2[[2]],1},{intersection2[[1]],intersection2[[2]],0}}]}];
point=Graphics3D[{PointSize->Large,Blue,Point[Append[a,0]]}];
tangent2=ParametricPlot[{t,a3D[[3]]+directionalDerivative[f,s,a3D]*t},
{t,-10,10},PlotStyle->{Thickness[0.003],Black}];
cut2=ParametricPlot[{t,f[a3D[[1]]+s[[1]]*t,a3D[[2]]+s[[2]]*t]},
{t,intersections[[1]],intersections[[2]]},PlotStyle->{Thickness[0.003],Black,Dashed}];
plane2=Graphics[{FaceForm[GrayLevel[0.5]],EdgeForm[],Opacity[0.5],
Polygon[{{intersections[[1]],0},{intersections[[1]],1},{intersections[[2]],1},{intersections[[2]],0}}]}];
point2=Graphics[{PointSize->Large,Blue,Point[{0,0}]}];
Labeled[GraphicsGrid[{{Show[{plotf,tangent,horizontalLine,cut,plane,point},
ViewPoint->{3,3,2}]},
{Show[{tangent2,cut2,plane2,point2},PlotRange->{{-2,2},{0,1.5}},AxesStyle->Blue,ImageSize->300,Axes->{True,False},Ticks->False]}},
Spacings->{0,5}],
Pane[Row[{Pane[Row[{"(OverVector[a]) = (",a3D[[1]],", ",a3D[[2]],")"}]],Spacer[20],Pane[Row[{"(OverVector[s]) =(",N[Round[s[[1]],0.1]],", ",N[Round[s[[2]],0.1]],")"}]],Spacer[20],Pane[Row[{"((∂f)/(∂(Overscript[s, ->])))((OverVector[a])) = ",N[Round[directionalDerivative[f,s,a3D],0.01]]}]]}],ImageMargins->20],Top]],{ϕ,0,2*π},ContinuousAction->False]

. Constrained extrema

The animation illustrates how the gradient of the function f(x,y)=((x^2)-(y^2)) (e^((-(x^2)-(y^2))/2)) and the gradient of the constraint (((x-1))^2)+(((y-1))^2)-(1/2)=0 look at the points of the constraint. The function f has two local extremum points subjected to the constraint g. Only at these points, the gradients of f and g are proportional. 
The gradient of f is displayed as a blue arrow, whereas the gradient of g is displayed as a red arrow. The critical points are represented by pink points. The length and the direction of the arrows correspond to the length and the direction of the gradient.

f[x_,y_]=((x^2)-(y^2))*(E^((-(x^2)-(y^2))/2));

g[x_,y_]= (((x-1))^2)+(((y-1))^2)-(1/2);

ax=1; ay=1; radius=(Sqrt[0.5]);

a=NSolve[g[x,y]==0&&D[f[x,y]-λ*g[x,y],x]==0 &&D[f[x,y]-λ*g[x,y],y]==0 ,{x,y,λ},Reals]

a1={x/.a[[1]],y/.a[[1]]};

a2={x/.a[[2]],y/.a[[2]]};

arrows[f_,startPoint_,colour_]:=
Module[{gradient,endPoint,scaleParameter},
scaleParameter=1;
gradient=Grad[f[x,y],{x,y}]/.{x->startPoint[[1]],y->startPoint[[2]]};
endPoint={gradient[[1]]*scaleParameter+startPoint[[1]],gradient[[2]]*scaleParameter+startPoint[[2]]};
Graphics[{colour,Arrow[{startPoint,endPoint}]}]]

backGround={ContourPlot[f[x,y],{x,-2,4},{y,-2,4},
Contours->10,ColorFunction->"BlueGreenYellow",
MaxRecursion->5,ImageSize->250,AxesLabel->{x,y}],
Graphics[Circle[{ax,ay},radius]],
Graphics[{PointSize[Large],Pink,Point[{a1[[1]],a1[[2]]}]}],
Graphics[{PointSize[Large],Pink,Point[{a2[[1]],a2[[2]]}]}]};

frames=Table[Show[Append[backGround,
{arrows[f,{Cos[angle]*radius+ax,Sin[angle]*radius+ay},Blue],arrows[g,{Cos[angle]*radius+ax,Sin[angle]*radius+ay},Red]}],
PlotLabel->Pane[Row[{"f[",N[Round[Cos[angle]*radius+ax,0.01]],", ",
N[Round[Sin[angle]*radius+ay,0.01]],"] = ",N[Round[f[Cos[angle]*radius+ax,Sin[angle]*radius+ay],0.01]]}]]],{angle,0,2*π,(π/50)}];

Labeled[ListAnimate[Rasterize[#,"Image"]&/@frames,AnimationRunning->False],PointLegend[{Pink,Blue,Red},{"critical points","gradient of f","gradient of the constaint"},LegendMarkers->{Graphics[Point[{0,0}]],"⟶" ,"⟶" }],Right]

Clear[f,g,ax,ay,radius,a,a1,a2,arrows,backGround,frames];

. Useful alternative coordinate systems

Transformation from cartesian to one of the following coordinate systems can help us with calculating integrals when integrated function and/or integration domain have radial, spherical or cylindrical symmetry.

.. Polar coordinate system

This is a two-dimensional coordinate system in which each point is determined by the distance from the origin and the angle from the polar axis. The relationship between these coordinate systems is given by

x=ρ cos(ϕ),
y=ρ sin(ϕ).

Panel[DynamicModule[{p={5,5}},Column[{Row[{Panel[Row[{"ρ = ",
Dynamic[EuclideanDistance[{0,0},p]]}],ImageSize->{120,50},Alignment->{Center,Center}],Panel[Row[{"ϕ = ",Dynamic[If[p[[2]]>0,
((VectorAngle[{1,0},p])/π),((2*π-VectorAngle[{1,0},p])/π)]]," π"}],
ImageSize->{120,50},Alignment->{Center,Center}]}],
Show[PolarPlot[{10},{t,0,2*π},PlotStyle->Lighter[Gray],PolarGridLines->{{0,Pi/2,Pi,3Pi/2},{2,4,6,8}}],Graphics[{Locator[
Dynamic[p]],Dynamic[Line[{{0,0},p}]],Opacity[0.2]},Axes->True],PlotRange->{{-11,11},{-11,11}},ImageSize->250]},Alignment->Center]]]

.. Spherical coordinate system

This coordinate system is three-dimensional. The position of a point is specified by the distance from the pole, the polar angle and the azimuth angle. The cartesian coordinates are determined by

x=ρ cos(θ) cos(ϕ),
y=ρ cos(θ) sin(ϕ),
z=ρ sin(θ).

Manipulate[
point={ρ*Cos[angles[[2]]]*Cos[angles[[1]]],ρ*Cos[angles[[2]]]*Sin[angles[[1]]],ρ*Sin[angles[[2]]]};
Show[Graphics3D[{PointSize[0.015],Point[point],
Thickness[0.003],Line[{{0,0,0},point}],Opacity[0],
EdgeForm[{Lighter[Gray],Thickness[0.001]}],
Cylinder[{{0,0,0},{0,0,0.01}},ρ],
EdgeForm[{Black}],
Cylinder[{{0,0,0},{((ρ*Cos[angles[[2]]]*Cos[angles[[1]]+(π/2)])/100),((ρ*Cos[angles[[2]]]*Sin[angles[[1]]+(π/2)])/100),0}},ρ],
EdgeForm[{Black}],
Cylinder[{{0,0,point[[3]]},{0,0,point[[3]]+0.01}},Abs[ρ*Cos[angles[[2]]]]],
Lighter[Gray],Opacity[0.15],Sphere[{0,0,0},ρ]},ViewPoint->{1,0.5,1}],PlotRange->{{-11,11},{-11,11},{-11,11}},ImageSize->{250},AxesOrigin->{0,0,0},Axes->True,Ticks->All,Boxed->False,BoxRatios->Automatic,SphericalRegion->True],
Column[{Row[{Control[{ρ,5,10}],Spacer[20],"ρ = ",Dynamic[ρ]}],
Row[{Labeled[Control[{{angles,{0,0},""},{0,-(π/2)},{2*π,(π/2)},ImageSize->Small}],{"ϕ","θ"},{Bottom,Left}],Spacer[160],Column[{Row[{"ϕ = ",Dynamic[((angles[[1]])/π)]," π"}],Row[{"θ = ",Dynamic[((angles[[2]])/π)]," π"}]}]}]}]]

Clear[point];

.. Cylindrical coordinate system

The cylindrical coordinate system is another three-dimensional coordinate system which is obtained by replacing x and y with polar coordinates ρ and ϕ and leaving the z coordinate the same, i.e.

x=ρ cos(ϕ)
y=ρ sin(ϕ)
z=h

Manipulate[
point={ρ*Cos[ϕ],ρ*Sin[ϕ],h};
Show[Graphics3D[
{PointSize[0.015],Point[point],
Thickness[0.003],Line[{{0,0,h},point}],
Lighter[Gray],Opacity[0.15],Cylinder[{{0,0,0},{0,0,h}},ρ]},
ViewPoint->{1,1,0.5}],PlotRange->{{-11,11},{-11,11},{-11,11}},ImageSize->{250},AxesOrigin->{0,0,0},Axes->True,Ticks->All,Boxed->False,BoxRatios->Automatic,SphericalRegion->True],
Column[{Row[{Control[{ρ,1,7}],Spacer[40],"ρ = ",Dynamic[ρ]}],
Row[{Control[{ϕ,0,2*π}],Spacer[40],"ϕ = ",Dynamic[(ϕ/π)]," π"}],
Row[{Control[{h,-7,7}],Spacer[40],"h = ",Dynamic[h]}]}]]

Clear[point];

. Double integral

This illustration shows how a definite integral of f(x,y)=2-2sin((x^2)+(y^2)) can be calculated using different numerical methods and the convergence of Riemann sums.
The interval of integration is divided into subintervals according to the number of dividing points. The approximation is performed by finding the volume of the collection of cuboids. The lengths and widths of these cuboids are defined by the size of subintervals. Their heights are determined by the approximation method:

The error is increasing or decreasing depending on the number of dividing points and the technique of approximation.

Manipulate[Module[{f,plotf,plotInterval,plotSize,intervalSize,partition,intervals,intervals2,fvalue,sum},
f[x_,y_]=2-2*Sin[(x^2)+(y^2)];
plotf=Plot3D[f[x,y],{x,-1,1},{y,-1,1},PlotStyle->{Opacity[0.5],Green},Mesh->None,Boxed->False,Axes->None];
plotInterval={-1,1};
plotSize=plotInterval[[2]]-plotInterval[[1]];
intervalSize=(plotSize/dividingPoints);
partition=Partition[Table[i,{i,plotInterval[[1]],plotInterval[[2]],intervalSize}],2,1];
intervals=Flatten[Table[{partition[[i]][[1]],partition[[j]][[2]]},{i,1,dividingPoints},{j,1,dividingPoints}],1];
intervals2=Flatten[Table[partition[[i]][[1]]<x<partition[[i]][[2]]&&partition[[j]][[1]]<y<partition[[j]][[2]],{i,1,dividingPoints},{j,1,dividingPoints}]];
fvalue=Switch[method,
lowerSum,Table[NMinimize[{f[x,y],intervals2[[i]]},{x,y},PrecisionGoal->2][[1]],{i,1,(dividingPoints^2)}],
upperSum,Table[NMaximize[{f[x,y],intervals2[[i]]},{x,y},PrecisionGoal->2][[1]],{i,1,(dividingPoints^2)}],
midpointRule,Table[f[intervals[[i]][[1]]+intervalSize/2,intervals[[i]][[2]]-intervalSize/2],{i,1,(dividingPoints^2)}],
trapezoidalRule,Table[
(1/4)*(f[intervals[[i]][[1]],intervals[[i]][[2]]]+f[intervals[[i]][[1]]+intervalSize,intervals[[i]][[2]]]+
f[intervals[[i]][[1]],intervals[[i]][[2]]-intervalSize]+f[intervals[[i]][[1]]+intervalSize,intervals[[i]][[2]]-intervalSize]),{i,1,(dividingPoints^2)}]];
sum=N[Total[Table[(intervalSize^2)*fvalue[[i]],{i,1,(dividingPoints^2)}]]];
Labeled[Show[{plotf,
Table[Graphics3D[{Opacity[0.7],Lighter[Yellow],Cuboid[{intervals[[i]][[1]],intervals[[i]][[2]],fvalue[[i]]},{intervals[[i]][[1]]+intervalSize,intervals[[i]][[2]]-intervalSize,0}]}],{i,1,(dividingPoints^2)}]},
ImageSize->{250,250},BoxRatios->{1,1,1},Boxed->False],
Pane[Column[{StandardForm[Panel[Row[{"f[x,y] = ",f[x,y]}]]],
StandardForm[Panel[Row[{"(∫_(-1)^1)(∫_(-1)^1) f[x,y] \[DifferentialD]x\[DifferentialD]y = ",N[Integrate[f[x,y],{x,-1,1},{y,-1,1}]]}]]],
Pane[Row[{"Volume = ",sum}],ImageMargins->20],
Pane[Row[{"Error = ",Abs[sum-N[Integrate[f[x,y],{x,-1,1},{y,-1,1}]]]}]]},
Center],ImageMargins->10],Right]],
{{dividingPoints,3},{3, 4,5,6,7,8,9,10}},{{method,lowerSum},{lowerSum, upperSum,midpointRule ,trapezoidalRule }},ControlType->{PopupMenu,PopupMenu}]

Bivariate limits

. Bivariate limits in Mathematica

Mathematica does not provide any function for calculating limits of functions of two variables. The only way is to use the command WolframAlpha which sends the input to WolframAlpha and imports the output. However, it requires the Internet access.

WolframAlpha["limit (x^2*y^2)/(x^2+y^2) as x->0, y->0","PodCells"]

Although it looks like it works pretty well, the results of some computations can be incorrect. For instance, the limit of the function ((2xy)/(xy + 2x - y)), where both x and y tend to zero.

WolframAlpha["limit (2*x*y)/(x*y + 2*x - y) as x->0, y->0","PodCells"]

WolframAlpha returns 0 but the limit does not exist as we proved in example ..

. Function doubleLimit

Due to the insufficient performance of the built-in function limit of WolframAlpha, we have decided to define our own function named doubleLimit. It consists of three modules, each of which uses a different method.
At first, we use the method of sequential limits. If both are finite, this module can determine that either the limit does not exist or the limit would equal some specific value if it existed.

sequentialLimits[f_,{x_,y_}->{ax_,ay_}]:=
Module[{seqLim2},
If[Head[f]===If||Head[f]===Piecewise||Head[f]===Which||Head[f]===Switch||Head[f]===Boole,"The limit may still exist.",
seqLim1=Limit[Limit[f,y->ay],x->ax];
seqLim2=Limit[Limit[f,x->ax],y->ay];
If[NumericQ[seqLim1]&&NumericQ[seqLim2]&&seqLim1≠seqLim2,"The limit does not exist.",
If[NumericQ[seqLim1]&&NumericQ[seqLim2]&&seqLim1==seqLim2,Row[{"The limit may exist and if it does then it has to be ",seqLim1,"."}],"The limit may still exist."]]]];

The second method how the nonexistence of a limit can be proven is approaching the point at which the limit is searched via various curves of the form x=(ky^(i/j)). If the function is a polynomial, the curves are produced for i and j run from 1 to the highest exponent from all of the numerators and denominators of the variables. In case the function is not polynomial, the highest exponent is increased by 1 in order to construct lines at least. The output of this module can be that the limit does not exist or what the limit would be provided it existed.

curves[f_,{x_,y_}->{ax_,ay_}]:=
Module[{g,nonexistence,exponents,highestExp,k},
If[Head[f]===If||Head[f]===Piecewise||Head[f]===Which||Head[f]===Switch||Head[f]===Boole,"The limit may still exist.",
g=Factor[f];
nonexistence=False;
exponents={Numerator[Exponent[Numerator[g],{x,y}]],Denominator[Exponent[Numerator[g],{x,y}]],
Numerator[Exponent[Denominator[g],{x,y}]],Denominator[Exponent[Denominator[g],{x,y}]]};
highestExp=Max[exponents]+1;
For[i=1,i≤highestExp,i++,
For[j=1,j≤highestExp,j++,
lim=Limit[g/.x->k*((Abs[y-ay])^(i/j))+ax,y->ay];
If[((Exponent[Numerator[lim],k]≠0||Exponent[Denominator[lim],k]≠0)&&Exponent[lim,k]≠-Infinity),{nonexistence=True,Break[]},Unevaluated[Sequence[]]]]];
If[nonexistence,"The limit does not exist.",If[NumericQ[lim],Row[{"The limit may exist and if it does then it has to be ",lim,"."}],"The limit may still exist."]]]];

The purpose of the last module is to find the limit employing the definition. Nonetheless, it can also state whether or not it exists. The function is maximized and minimized on two-dimensional interval of the size 2d × 2d, where the point at which the limit should be found is in the middle of this interval. The result of this maximization and minimization are one-variable functions dependent on d. If they have the same limit as d tends to zero, the searched limit exists and is equal to the their common value. If they differ, the limit does not exist. This method, which is described at [], works fine typically for rational functions.

maxMinMethod[f_,{x_,y_}->{ax_,ay_}]:=
Module[{d,max,min,limMax,limMin},
max=Maximize[{f,0<d<1&&-d≤x-ax≤d&&-d≤y-ay≤d},{x,y}];
min=Minimize[{f,0<d<1&&-d≤x-ax≤d&&-d≤y-ay≤d},{x,y}];
limMax=Limit[max[[1]],d->0];
limMin=Limit[min[[1]],d->0];
If[(limMax==Infinity && limMin==-Infinity)||(limMax==-Infinity && limMin==Infinity),"The limit does not exist.",
If[!NumericQ[limMax]||!NumericQ[limMin],"Function cannot find the limit.",
If[limMax==limMin,limMax,"The limit does not exist."]]]];

Note that this method can be easily transformed to search limits of more than two-variable functions.
All three methods are used in the following complex module termed doubleLimit.

doubleLimit[f_,{x_,y_}->{ax_,ay_}]:=
Module[{possibleLimit,limit,seqLim1,lim},
possibleLimit=False;
limit=sequentialLimits[f,{x,y}->{ax,ay}];
If[limit==="The limit does not exist.",limit,If[limit===Row[{"The limit may exist and if it does then it has to be ",seqLim1,"."}],possibleLimit=limit];
limit=curves[f,{x,y}->{ax,ay}];
If[limit==="The limit does not exist.",limit,If[limit===Row[{"The limit may exist and if it does then it has to be ",lim,"."}]
&&possibleLimit===False,possibleLimit=limit];
limit=maxMinMethod[f,{x,y}->{ax,ay}];
If[limit==="The limit does not exist."||NumericQ[limit],limit,If[possibleLimit=!=False,possibleLimit,"Function cannot find the limit."]]]]];

Here are some examples of applications:

doubleLimit[(((x^2)*(y^2))/((x^2)+(y^2))),{x,y}->{0,0}]

doubleLimit[((2*x*y)/(x*y+2*x-y)),{x,y}->{0,0}]

doubleLimit[(((x^2)-(y^2))/((x^2)+(y^2))),{x,y}->{0,0}]

doubleLimit[Cos[π*((x^2)/(4*(x^2)+(y^4)))],{x,y}->{0,0}]

Clear[sequentialLimits,seqLim1,curves,lim,maxMinMethod,doubleLimit];

.. Method of critical paths

There is one more method published at [] which we intended to employ for calculating bivariate limits. It is used and delineated by Maplesoft in Maple version 17. It is based on so called critical paths. The function attains its local maximum or minimum subjected to the constraint C, which is a circle with radius ρ and centre (OverVector[a])=((Subscript[a, x]), (Subscript[a, y])), if the gradient of the function and the gradient of the constraint are parallel. Thus, the following equation can be derived:

((∂f)/(∂x))((∂C)/(∂y))-((∂f)/(∂y))((∂C)/(∂x))=0.

Together with the equation of the function, it determines the critical paths. If all limits along these paths are the same, the bivariate limit exists and equals this common value. On the other hand, if at least one of them differs, it can be said that the bivariate limit does not exist.
Although it works in theory, the reality is much worse. Function Solve has only limited power how to find the solution, especially of more complicated equations. Moreover, this method solves the same kind of examples as the module maxMinMethod.
Taking into consideration these facts, we have decided not to use this method.

.. Possible issues

The success of finding the bivariate limit depends on the built-in function Limit which does not have to be so trustworthy as one would expect. For example, the following limit should not exist as the one-sided  limits differ.

Limit[((Abs[x])/x),x->0]

However, Mathematica declares that it exists because, by default, it computes the limit from the right.

Ballistic curves

A ballistic curve is a path of a projectile in the gravity field. Besides the gravitational force, there are other forces acting on the projectile. One of the them is air resistance (also called drag). The earth's rotation is neglected, and the height between the projectile and the earth is assumed to be insignificant in comparison with the earth's dimensions in order that the gravitational constant would not change.

. Ballistic trajectory without air resistance

The projectile is launched at a velocity of

(OverVector[v])=((Subscript[(Subscript[v, x]), ,])(Subscript[v, y])) .

The x- and y-component can be expressed as

(Subscript[v, x])=v cos(ϕ) ,

(Subscript[v, y])=v sin(ϕ),

where ϕ is the angle between x-axis and vector (OverVector[v]).
Let us assume that the projectile is launched from the point (0, 0). The gravitational force influences only the vertical component. The trajectory is given by the parametric equations

x(t)=(Subscript[v, x]) t,

y(t)=(Subscript[v, y]) t-(1/2)(gt^2).

. Ballistic trajectory with air resistance

Air resistance is a kind of force which acts in the opposite way to the motion of an object with respect to the surrounding fluid [].
We will consider two kinds of drag ‒ linear and quadratic. The former is usually used for lower velocity, the latter usually for higher velocity []. They are given by the equations

(Subscript[(OverVector[F]), L])=-c(OverVector[v]), []

(Subscript[(OverVector[F]), Q])=-c ((OverVector[v])^2), []

where c is the drag coefficient that is used to quantify the drag of an object. It is associated with the surface area of the object and is obtained from laboratory experiments. The value usually ranges between 0.4 and 1 [].

.. Drag proportional to velocity

The total force in x- and y-direction should be expressed as

(Subscript[F, x])=(Subscript[ma, x])=-(Subscript[cv, x])

and

(Subscript[F, y])=(Subscript[ma, y])=-(Subscript[cv, y])-mg.

Note that only (Subscript[F, y]) is influenced by the gravitational force. Since the trajectory should be found, the following differential equations are derived by alteration of () and ():

(((d^2)x)/(d(t^2)))=-(c/m)((dx)/(dt)),

(((d^2)y)/(d(t^2)))=-(c/m)((dy)/(dt))-g.

As equations () and () are independent, they can be solved separately. The particular solutions are fixed by the initial conditions x'(0)=(Subscript[v, x]), y'(0)=(Subscript[v, y]), x(0) = 0, and y(0)=0.

DSolve[{x''[t]==-c*((x'[t])/m),x'[0]==vx,x[0]==0},x[t],t]//FullSimplify

DSolve[{y''[t]==-c*((y'[t])/m)-g,y'[0]==vy,y[0]==0},y[t],t]//FullSimplify

The equations of the trajectory are

x(t)=(((1-(e^(-((c t)/m)))) m (Subscript[v, x]))/c),

y(t)=((m (g m-c g t+c (Subscript[v, y])-(e^(-((c t)/m))) (g m+c (Subscript[v, y]))))/(c^2)).

.. Drag proportional to velocity squared

The total force in x- and y-direction is

(Subscript[F, x])=(Subscript[ma, x])=-((Subscript[cv, x])^2),

(Subscript[F, y])=(Subscript[ma, y])=-((Subscript[cv, y])^2)-mg.

Thus, the trajectory can be expressed as the following differential equations:

(((d^2)x)/(d(t^2)))=-(c/m)((dx)/(dt))(Sqrt[(((((dx)/(dt))))^2)+(((((dy)/(dt))))^2)]),

(((d^2)y)/(d(t^2)))=-(c/m)((dy)/(dt))(Sqrt[(((((dx)/(dt))))^2)+(((((dy)/(dt))))^2)])-g.

Unlike () and (), equations () and () are dependent and must be solved as a system of differential equations. However, Mathematica is not able to solve them analytically. As a consequence, NDSolve is used instead of DSolve to find numerical solution. Hence, c, m, and g must be equal to some specific value. The initial conditions are x'(0)=(Subscript[v, x]), y'(0)=(Subscript[v, y]), x(0) = 0, and y(0)=0.

NDSolve[{x''[t]==-(c/m)*x'[t]*(Sqrt[((x'[t])^2)+((y'[t])^2)])/.{c->1,m->1,g->9.81},y''[t]==-g-(c/m)*y'[t]*(Sqrt[((x'[t])^2)+((y'[t])^2)])/.{c->1,m->1,g->9.81},x'[0]==1,y'[0]==1,x[0]==0,y[0]==0},{x[t],y[t]},{t,0,3}]

NDSolve returns only an approximate function, which differs according to the specific initial conditions.

. Visualization

The trajectory of a projectile is now visualized using the equations () and () for trajectory without drag, and () and () for trajectory with drag proportional to the velocity. Drag proportional to the velocity squared is determined numerically according to the selected initial conditions. The mass m and constant c must be fixed for the purposes of visualization (we chose m=1 and c=1).

Manipulate[Module[{x1,y1,x2,y2,numSol},
x1[t_]=vx*t;
y1[t_]=vy*t-(1/2)*g*(t^2);
If[f1==1,plot1=ParametricPlot[
{{x1[t]/.{vx->v0*Cos[ϕ0],ϕ->ϕ0,m->1,g->9.81,c->1},y1[t]/.{vy->v0*Sin[ϕ0],ϕ->ϕ0,m->1,g->9.81,c->1}}},{t,0,3},PlotStyle->Blue],plot1=ParametricPlot[0,{t,0,3}]];
x2[t_]=(((1-(E^(-((c*t)/m))))*m*vx)/c);
y2[t_]=((m*(g*m-c*g*t+c*vy-(E^(-((c*t)/m)))*(g*m+c*vy)))/(c^2));
If[f2==1,plot2=ParametricPlot[
{{x2[t]/.{vx->v0*Cos[ϕ0],ϕ->ϕ0,m->1,g->9.81,c->1},y2[t]/.{vy->v0*Sin[ϕ0],ϕ->ϕ0,m->1,g->9.81,c->1}}},{t,0,3},PlotStyle->Red],plot2=ParametricPlot[0,{t,0,3}]];
numSol=NDSolve[{x3''[t]==-(c/m)*x3'[t]*(Sqrt[((x3'[t])^2)+((y3'[t])^2)])/.{c->1,m->1,g->9.81},y3''[t]==-g-(c/m)*y3'[t]*(Sqrt[((x3'[t])^2)+((y3'[t])^2)])/.{c->1,m->1,g->9.81},x3[0]==0,y3[0]==0,x3'[0]==v0*Cos[ϕ0],y3'[0]==v0*Sin[ϕ0]},{x3[t],y3[t]},{t,0,3}];
If[f3==1,plot3=ParametricPlot[{x3[t],y3[t]}/.numSol,{t,0,3},PlotStyle->Darker[Green]],plot3=ParametricPlot[0,{t,0,3}]];
Show[{plot1,plot2,plot3},PlotRange->{{0,11},{0.1,6}},ImageSize->400]],
Column[{Control[{{v0,7.5,"v"},5,10}],Control[{{ϕ0,(π/4),"ϕ"},0,(π/2)}],
Row[{Control[{{f1,1,""},{0,1}}],Text[Style[ " without drag",Medium,Blue]]}],
Row[{Control[{{f2,1,""},{0,1}}],Text[Style[ " drag proportional to velocity",Medium,Red]]}],
Row[{Control[{{f3,1,""},{0,1}}],Text[Style[ " drag proportional to velocity squared",Medium,Darker[Green]]]}]}]]

Conclusion

An interactive tool for students of the Mathematics 4 course has been developed. This project could help them to understand the problems on many solved examples. The advantages of computer algebra system Wolfram Mathematica have been demonstrated on particular examples. On the other hand, we were also able to show that Mathematica is not unfailing and sometimes gives bad results without any warning.
Since the source code is included for all examples, this project can also help learning the basics of Mathematica, although it is not the main purpose.
The question of bivariate limits has been discussed. Our own function has been created to find them or disprove their existence using various methods. It works very good; especially when it is applied to a rational function, and covers the needs of the Mathematics 4 course.
Furthermore, the problem of ballistic curves has been solved and visualized.
Several interactive visualizations have been developed to help students understand some notions and problems. Moreover, they could be also used as an interactive tool by lecturers.
In the future, the function for calculating bivariate limits could be upgraded to operate on more difficult examples. Furthermore, other courses of mathematics could be processed in a similar way to create complex materials for mathematics courses at our Faculty.

Bibliography

[]		Abell, Martha L., and James P. Braselton. Mathematica By Example. 3rd ed. Boston: Academic, 2004. Print.

[]		Halliday, David, Robert Resnick, and Jearl Walker. Fyzika: Část 1 Mechanika. Trans. Petr Dub and Jan Obdržálek. Brno: VUTIUM, 2000. Print.

[]		Mangano, Sal. Mathematica Cookbook. Sebastopol, CA: O'Reilly, 2010. Print.

[]		Marsden, Jerrold and Alan Weinstein. Calculus III. New York: Springer-Verlag, 1985. Print.

[]		Štoll, Ivan. Mechanika. Prague: Česká Technika, 2010. Print.

[]		About Wolfram|Alpha. Making the World's Knowledge Computable. N.p., [2014]. Web. 24 June 2014. <http://www.wolframalpha.com/about.html>.

[]		About Wolfram Research. Wolfram Research Company Background. N.p., [2014]. Web. 24 June 2014. <http://www.wolfram.com/company/background.html>.

[]		Bivariate Limits. Updates/Maple17/BivariateLimits - Help. N.p., [2014]. Web. 25 June 2014. <http://www.maplesoft.com/support/help/Maple/view.aspx?path=updates/Maple17/BivariateLimits>.

[]		Dawkins, Paul. Pauls Online Notes : Calculus II - Approximating Definite Integrals. N.p., [2014]. Web. 24 June 2014. <http://tutorial.math.lamar.edu/Classes/CalcII/ApproximatingDefIntegrals.aspx>.

[]		Erlichson, Herman. Maximum Projectile Range with Drag and Lift, with Particular Application to Golf (n.d.): n. pag. 16 June 1982. Web. 24 June 2014. <http://ww2.odu.edu/~agodunov/teaching/phys420/files/Erlichson.pdf>.

[]		Franc, Jiří, and Matěj Tušek. Sbírka úloh k předmětu Matematika 4 (MAT4) (n.d.): n. pag. 25 Feb. 2014. Web. 24 June 2014. <http://kmlinux.fjfi.cvut.cz/~tusekmat/download/MAT4_sbirka.pdf>.

[]		Michael E2. Finding Limits in Several Variables. Functions. N.p., 18 Mar. 2013. Web. 24 June 2014. <http://mathematica.stackexchange.com/questions/21544/finding-limits-in-several-variables>.

[]		Wikipedia Contributors. Drag (physics). Wikipedia. Wikimedia Foundation, 18 June 2014. Web. 24 June 2014. <http://en.wikipedia.org/wiki/Drag_%28physics%29>.

[]		Wikipedia Contributors. Mathematica. Wikipedia. Wikimedia Foundation, 22 June 2014. Web. 25 June 2014. <https://en.wikipedia.org/wiki/Mathematica>.

[]		Wolfram|Alpha. Personal Analytics. N.p., [2014]. Web. 24 June 2014. <http://www.wolframalpha.com/facebook/>.

[]		Wolfram Demonstrations Project. Demonstrations RSS. N.p., [2014]. Web. 24 June 2014. <http://demonstrations.wolfram.com/>.