(* ::Package:: *) (* ::Input::Initialization:: *) (* ::Input::Initialization:: *) (*Usage*) VecPlot::usage="Plots a list of vectors on an i,j plane."; DetailBase::usage="Returns the components required for DetailPlot."; DetailPlot::usage = "Plots asymptotes and other important details of a graph."; RelationPlot::usage = "Plots relations returning axial intercepts."; Va::usage="Returns a vertical asymptote at x=a."; Euler::usage = "Runs Euler's method."; FindRoots::usage ="FindRoot but it attempts to find more than one root. Input is FindRoots[function,{variable,min,max}]" ImpDiff::usage = "Implicitly differentiates an expression once or twice. Input is ImpDiff[expression,{x,y,number of difs}]" TangentLine::usage=" TangentLine[function,variable,point] returns the equation of the tangent at a point." ForceSolve::usage = "Forces Mathematica to return a numerical expression from a root expression."; SumToProduct::usage = "Converts the product of sine and cosine functions to the sum of cosine and sine functions. Usage is SumToProduct[function1,function2] if your function is function1*function2."; ToDeg::usage = "Converts radians to degrees."; ToRad::usage="Converts degrees to radians."; Compsq::usage="Convert standard form quadratic to turning point form."; HGD::usage="Shortcut for HypergeometricDistribution. Input is number of trials, defectives and total population"; ND::usage="Shortcut for NormalDistribution. Input is \[Mu] and \[Sigma]."; HeronArea::usage="Returns the area of a triangle when three sides, a,b and c are known, Input is HeronArea[a,b,c]"; Dp::usage = "Gives a number to n decimal places. Input is Dp[expression,n]"; ScalarRes::usage = "Gives the scalar resolute of u onto v. Usage is ScalarRes[u,v]"; StrictFunctionRange::usage "\ StrictFunctionRange[expr, var1,var2] \ finds the largest range of definition of expression expr, treated as function \ in given variables vars. \ Excludes all singularities encountered while evaluating expr, \ even if those singularities would be removed \ during ordinary evaluation of expr. \ var1 and var2 can be a symbol or list of symbols. Domain dom can be Reals or Complexes."; StrictFunctionDomain::usage = "\ StrictFunctionDomain[expr, vars, dom] \ finds the largest domain of definition of expression expr, treated as function \ in given variables vars, with arguments and values in domain dom, \ excludes all singularities encountered while evaluating expr, \ even if those singularities would be removed \ during ordinary evaluation of expr. \ vars can be a symbol or list of symbols. Domain dom can be Reals or Complexes.\ StrictFunctionDomain[expr, vars] \ uses Reals as domain.\ StrictFunctionDomain[expr] or StrictFunctionDomain[expr, Automatic, ...] \ uses variables extracted from expr."; RestrictDomain::usage = "\ RestrictDomain[expr, vars, dom] \ returns expr with sub-expressions subExpr, for which, \ domain can be restricted, replaced with ConditionalExpression[subExpr, cond] \ where cond are conditions restricting domain of subExpr, \ treated as function in given variables vars, \ with arguments and values in domain dom. \ vars can be a symbol or list of symbols. Domain dom can be Reals or Complexes.\ RestrictDomain[expr, vars] \ uses Reals as domain.\ RestrictDomain[expr] or RestrictDomain[expr, Automatic, ...] \ for each replaced sub-expression uses variables extracted from it."; Begin["`Private`"]; SetOptions[$FrontEndSession, InputAutoReplacements ->{ "m11"->ToBoxes[Defer[FullSimplify/@{}]], "m12"->ToBoxes[Defer[Solve[a,x]]], "m13"->ToBoxes[Defer[Reduce/@{}]], "m14"->ToBoxes[Defer[ComplexExpand/@{}]], "m15"->ToBoxes[Defer[ArcLength[{x[t],y[t]},{t,a,b}]]], "m16"->ToBoxes[Defer[Euler[exp,{x,y},{Subscript[x, o],Subscript[y, o]},h,10]]], "m17"->ToBoxes@Defer[Solve[D[y,x],y'[x]]], "m18"->ToBoxes[Defer[Det[{{},{},{}}]]], "m19"->ToBoxes[Defer[Integrate[y,{x,a,b}]]], "m21"->ToBoxes[Defer[Plot[fn,{x,-10,10},PlotRange->10]]], "m22"->ToBoxes[Defer[ContourPlot[rel,{x,-10,10},{y,-10,10},Frame->False,Axes->True]]], "m23"->ToBoxes[Defer[StreamPlot[{1,dy/dx},{x,-10,10},{y,-10,10},Frame->False,Axes->True]]], "m23a"->ToBoxes[Defer[StreamPlot[{1,#},{x,-10,10},{y,-10,10},Frame->False,Axes->True]&/@{}]], "m24"->ToBoxes[Defer[ParametricPlot[{x[t],y[t]},{t,a,b}]]], "m25"->ToBoxes[Defer[DetailPlot[fn,{x,-10,10},PlotRange->10]]], "m25a"->ToBoxes[Defer[DetailPlot[#,{x,-10,10},PlotRange->10]&/@{}]], "m26"->ToBoxes[Defer[PolarPlot[r[\[Theta]],{\[Theta],0,2\[Pi]}]]], "m27"->ToBoxes[Defer[RelationPlot[rel,{x,-5,5},{y,-5,5}]]], "m27a"->ToBoxes[Defer[(RelationPlot[#,{x,-5,5},{y,-5,5}])&/@{}]] }] (*Graphing Utilities*) Va[x_]:=InfiniteLine[{x,0},{0,-x}] Options[DetailBase]=Options[Plot]; DetailBase[exp_,{x_,min_,max_},opts:OptionsPattern[]]:= Module[{xint,xintP,yint,yintP,stat,statP,expD,inflec,inflecP, term,residue,obliq,vert,line,vaP,vaN,obliqStr,plt,domain,holex,holeP,va,vaC,oblf,asym,cond}, expD=D[exp,x]; (*Points*) term=If[Length[exp]>1,Level[Apart@exp,1],Apart@exp]; domain=FunctionDomain[{exp,min<=x<=max},x]; xint=Thread[{ForceSolve[{exp==0,domain},x],0}];(*limit domain for trig functions*) xintP=Graphics[{PointSize[0.015],Red,#}]&/@(Tooltip[Point[#],#]&/@xint); cond=If[FunctionDomain[exp,x]==True,True,False]; holex=Complement[Thread[{ForceSolve[{TrigReduce@exp==0,min<=x<=max},x],0}],xint]; holeP=Graphics[{PointSize[0.015],White,#}]&/@(Tooltip[Point[#],#]&/@holex); yint=Solve[y==exp&&x==0&&domain,{x,y},Reals][[All,All,2]]; yintP=Graphics[{PointSize[0.015],Red,#}]&/@(Tooltip[Point[#],#]&/@yint); stat=Thread[{ForceSolve[{expD==0,domain},x],ReplaceAll[exp,x->#]&/@ForceSolve[{expD==0,domain},x]}];statP=Graphics[{PointSize[0.015],Blue,#}]&/@(Tooltip[Point[#],#]&/@stat);inflec=Thread[{ForceSolve[{D[expD,x]==0,domain},x],ReplaceAll[exp,x->#]&/@ForceSolve[{D[expD,x]==0,domain},x]}];inflecP=Graphics[{PointSize[0.015],Green,#}]&/@(Tooltip[Point[#],#]&/@inflec); (*Asymptotes*)If[NumberQ@Limit[exp,x->\[Infinity]],vaP=Limit[exp,x->\[Infinity]],Nothing]; If[NumberQ@Limit[exp,x->-\[Infinity]],vaN=Limit[exp,x->-\[Infinity]],Nothing]; term=If[Length[exp]>1,Level[Apart@exp,1],Apart@exp]; residue=If[Length[exp]>1,Plus@@Select[term,(0==Limit[#,x->\[Infinity]]&&0==Limit[#,x->-\[Infinity]])&],exp]; obliq=If[!NumericQ[Limit[exp,x->\[Infinity]]],If[NumericQ[Apart@exp-residue],Nothing,Apart@exp-residue],Limit[exp,x->\[Infinity]]]; If[cond,obliq=Nothing]; (*Graphics*) vaC=Quiet@ContourPlot[Evaluate[Labeled/@ (*Here it takes only asymptotes that aren't also holes*)((Union[Complement[Flatten@DeleteCases[Solve[TrigReduce[Together[1/#]]==0&&min<=x<=max,x]&/@Flatten[{term}],{},Infinity],Flatten@Solve[{TrigReduce@exp==0,min<=x<=max},x]],Flatten@Solve[Denominator[exp]==0&&min<=x<=max,x]])/.Rule->Equal)],(*Then tries to solve for the reciprocal being 0*) {x,-10^4,10^4},{y,-10^4,10^4},ContourStyle->{Directive[Red,Dashed]}]; asym=Plot[Tooltip[{obliq,vaP,vaN}],{x,min,max},PlotStyle->Directive[Red,Dashed],Evaluate[FilterRules[{opts},{Except[PlotStyle],Options[Plot]}]]]; {asym,xintP,yintP,holeP,statP,inflecP,vaC}]; Options[DetailPlot]=Options[Plot]; DetailPlot[exp_,{x_,min_,max_},opts:OptionsPattern[]]:= Quiet@Module[{bases,asym,base,reVar}, If[ListQ[exp], bases=DetailBase[#,{x,min,max}]&/@exp; base=Plot[Tooltip[exp], {x,min,max}, Evaluate[FilterRules[{opts},Options[Plot]]]]; Show[Flatten[{base,bases}]], bases=DetailBase[exp,{x,min,max}]; base=Plot[Tooltip[exp], {x,min,max}, Evaluate[FilterRules[{opts},Options[Plot]]]]; Show[Flatten[{base,bases}]] ] ] Options[RelationPlot]=Options[ContourPlot]; RelationPlot[exp_,{x_,min_,max_},{y_,yin_,yang_},opts:OptionsPattern[]]:= Module[{ex,exf,xint,xintP,yint,yintP,stat,statP,expD,base,vert,vertP,diff,hor,horP}, exf=ComplexExpand[exp]; (*Points*) xint=Solve[exf&&y==0&&min<=x<=max,{x,y},Reals][[All,All,2]];(*limit domain for trig functions*) xintP=Graphics[{PointSize[0.015],Red,#}]&/@(Tooltip[Point[#],ToString[#,TraditionalForm]]&/@xint); yint=Solve[exf&&x==0,{x,y},Reals][[All,All,2]]; yintP=Graphics[{PointSize[0.015],Red,#}]&/@(Tooltip[Point[#],ToString[#,TraditionalForm]]&/@yint); diff=y'[x]/.Solve[D[exp/.y->y[x],x],y'[x]][[1]]/.y[x]->y; vert=Solve[Denominator[diff]==0&&exp&&min<=x<=max&&yin<=y<=yang,{x,y}][[All,All,2]]//ToRadicals; vertP=Graphics[{PointSize[0.015],Blue,#}]&/@(Tooltip[Point[#],ToString[#,TraditionalForm]]&/@vert); hor=Solve[Numerator[diff]==0&&exp&&min<=x<=max&&yin<=y<=yang,{x,y}][[All,All,2]]//ToRadicals; horP=Graphics[{PointSize[0.015],Blue,#}]&/@(Tooltip[Point[#],ToString[#,TraditionalForm]]&/@hor); (*Graphics*)base=ContourPlot[exp,{x,min,max},{y,yin,yang},Frame->False,Axes->True,Evaluate[FilterRules[{opts},Options[ContourPlot]]]]; Show[Flatten[{base,xintP,yintP,vertP,horP}]]] VecPlot[vecs__?ListQ]:= Module[{vec}, vec=Arrow[{{0,0},#}]&/@vecs; Graphics[vec,Axes->True, AxesLabel->{"i","j"},ImageSize->Large,Background->White]] (*Algebra*) TangentLine[f_,var_,a_]:=Expand[y/.Solve[y-(f/.var->a)==(D[f,var]/.var->a)(var-a),y]]; ScalarRes[u_,v_]:= Simplify@Dot[u,Normalize[v]]; ForceSolve[stuff_,variables_]:= Module[{b,c}, a=Solve[stuff,variables]; If[a=={{}},Return[{}]]; a=a[[All,1,2]]; b=DeleteCases[a,_Root]; c=N/@Cases[a,_Root]; Join[b,c] ]; (*Spreadsheet and iterative functions*) Euler[f_, {varx_, vary_}, {xo_, yo_}, h_, xF_] := TableForm@ Join[{{ToString[varx], ToString[vary]}}, NestWhileList[ Function[pointxy, pointxy + {h, h*Apply[Function[{varx, vary}, f], pointxy]}], {xo, yo}, If[h > 0, #[[1]] < xF &, #[[1]] > xF &]]]; Options@FindRoots=Sort@Join[Options@FindRoot,{MaxRecursion->Automatic,PerformanceGoal:>$PerformanceGoal,PlotPoints->Automatic,Debug->False,ZeroTolerance->10^-2}]; FindRoots[fun_,{var_,min_,max_},opts:OptionsPattern[]]:=Module[{PlotRules,RootRules,g,g2,pts,pts2,lpts,F,sol}, (*Extract the Options*)PlotRules=Sequence@@FilterRules[Join[{opts},Options@FindRoots],Options@Plot]; RootRules=Sequence@@FilterRules[Join[{opts},Options@FindRoots],Options@FindRoot]; (*Plot the function and "mesh" the point with y-coordinate 0*)g=Normal@Plot[fun,{var,min,max},MeshFunctions->(#2&),Mesh->{{0}},Method->Automatic,Evaluate@PlotRules]; (*Get the meshes zeros*)pts=Cases[g,Point[p_]:>SetPrecision[p[[1]],OptionValue@WorkingPrecision],Infinity]; (*Get all plot points*)lpts=Join@@Cases[g,Line[p_]:>SetPrecision[p,OptionValue@WorkingPrecision],Infinity]; (*Derive the interpolated data to find other zeros*)F=Interpolation[lpts,InterpolationOrder->2]; g2=Normal@Plot[Evaluate@D[F@var,var],{var,min,max},MeshFunctions->(#2&),Mesh->{{0}},Method->Automatic,Evaluate@PlotRules]; (*Get the meshes zeros and retain only small ones*)pts2=Cases[g2,Point[p_]:>SetPrecision[p[[1]],OptionValue@WorkingPrecision],Infinity]; pts2=Select[pts2,Abs[F@#]0,pts=Map[FindRoot[fun,{var,#},Evaluate@RootRules]&,pts]; sol=Union@Select[pts,min<=Last@Last@#<=max&]; (*For debug purposes*)If[OptionValue@Debug,Print@Show[g,Graphics@{PointSize@0.02,Red,Point[{var,fun}/.sol]}]]; sol,If[OptionValue@Debug,Print@g]; {}]] (*Function domain and range*) StrictFunctionRange // Attributes = HoldFirst; StrictFunctionRange[ expr_,var1_,var2_,dom:Real | Complexes : Reals, opts : OptionsPattern@FunctionRange ]:= With[{eval = evaluateKeepingSingularities@expr}, a = FunctionDomain @@ Unevaluated /@ Join[ eval, getHeldVars[var1]@eval, HoldComplete[dom, opts] ]; Reduce[{var2==expr,a},var2,{var1},Reals] ] StrictFunctionDomain // Attributes = HoldFirst; StrictFunctionDomain[ expr_, vars_ : Automatic, dom : Reals | Complexes : Reals, opts : OptionsPattern@FunctionDomain ] := With[{eval = evaluateKeepingSingularities@expr}, FunctionDomain @@ Unevaluated /@ Join[ eval, getHeldVars[vars]@eval, HoldComplete[dom, opts] ] ] RestrictDomain // Attributes = HoldFirst; RestrictDomain[ expr_, vars_ : Automatic, dom : Reals | Complexes : Reals, opts : OptionsPattern@FunctionDomain ] := ReleaseHold@With[{getHeldVars = getHeldVars@vars}, evaluateKeepingSingularities@expr /. subExpr : (acceptableHeads@dom)[___] :> ConditionalExpression[subExpr, FunctionDomain @@ Unevaluated /@ Join[ HoldComplete@subExpr, getHeldVars@subExpr, HoldComplete[dom, opts] ] ] ] evaluateKeepingSingularities = Function[Null, Internal`InheritedBlock[{Power}, With[{protected = Unprotect@Power}, Power@args__ /; Not@VectorQ[{args}, NumericQ] := power@args; Protect@protected ]; HoldComplete@#&@# /. power -> Power ], HoldAllComplete ]; getHeldVars // Attributes = HoldAllComplete; getHeldVars@Automatic = Function[Null, (HoldComplete[#]&@Union@Cases[Unevaluated@#, s_Symbol /; Not@MemberQ[Attributes@s, Constant] :> HoldComplete@s, {-1} ])[[All, All, 1]], HoldAllComplete ]; getHeldVars@vars_ := Function[Null, HoldComplete@vars, HoldAllComplete] acceptableHeads@Complexes = Abs | AiryAi | AiryAiPrime | AiryBi | AiryBiPrime | AngerJ | ArcCos | ArcCosh | ArcCot | ArcCoth | ArcCsc | ArcCsch | ArcSec | ArcSech | ArcSin | ArcSinh | ArcTan | ArcTanh | Arg | BesselI | BesselJ | BesselK | BesselY | Beta | Binomial | Boole | Ceiling | ConditionalExpression | Conjugate | Cos | Cosh | CoshIntegral | CosIntegral | Cot | Coth | Csc | Csch | CubeRoot | DawsonF | DiscreteDelta | Erf | Erfc | Erfi | Exp | ExpIntegralE | ExpIntegralEi | Factorial | Factorial2 | Fibonacci | Floor | FractionalPart | FresnelC | FresnelS | Function | Gamma | GammaRegularized | GegenbauerC | Gudermannian | HarmonicNumber | Haversine | HermiteH | Hypergeometric0F1 | Hypergeometric0F1Regularized | Hypergeometric1F1 | Hypergeometric1F1Regularized | Im | IntegerPart | Integrate | KroneckerDelta | LambertW | List | Log | Log10 | Log2 | LogGamma | LogIntegral | Max | Min | Mod | Norm | Plus | Pochhammer | PolyGamma | Power | PrimePi | ProductLog | Quotient | Re | RiemannSiegelTheta | RiemannSiegelZ | Round | SawtoothWave | Sec | Sech | Sign | Sin | Sinc | Sinh | SinhIntegral | SinIntegral | Sqrt | SquareWave | Surd | Tan | Tanh | Times | TriangleWave | UnitStep | WeberE | Zeta; acceptableHeads@Reals = Union[acceptableHeads@Complexes, BarnesG | EllipticE | EllipticK | EllipticNomeQ | InverseEllipticNomeQ | InverseErf | InverseErfc | InverseGammaRegularized ]; ImpDiff[expr_,{x_:x,y_:y[x],n_:1}]:= Solve[D[expr,{x,n}],D[y,{x,n}]] (* ::Input::Initialization:: *) (*quick function use*) SumToProduct[a_,b_]:= Switch[a[[0]], Cos,Switch[b[[0]], Sin,1/2 (Sin[a[[1]]+b[[1]]]-Sin[a[[1]]-b[[1]]]), Cos,1/2 (Cos[a[[1]]-b[[1]]]+Cos[a[[1]]+b[[1]]])], Sin,Switch[b[[0]], Sin,1/2 (Cos[a[[1]]-b[[1]]]-Cos[a[[1]]+b[[1]]]), Cos,1/2 (Sin[a[[1]]+b[[1]]]+Sin[a[[1]]-b[[1]]])]]; ToDeg[Radian_]:= N@Radian*180/Pi; ToRad[Degree_]:= N@Degree*Pi/180; Compsq[a_,b_,c_]:=a (x+b/(2a))^2+((4a*c-b^2)/(4a))//TraditionalForm; HGD[n_,D_,N_]:=HypergeometricDistribution[n,D,N]; ND[\[Mu]_,\[Sigma]_]:=NormalDistribution[\[Mu],\[Sigma]]; HeronArea[a_,b_,c_]:= Module[{s}, s = (a+b+c)/2; Sqrt[s(s-a)(s-b)(s-c)]]; Dp[expr_,n_]:= N@expr~NumberForm~n (* ::Input::Initialization:: *) End[] EndPackage[]