Re: Plotting two surfaces
- Subject: Re: Plotting two surfaces
- From: uunet!garnet.berkeley.edu!campbell (Robert I. Campbell)
- Date: Mon, 12 Mar 90 16:16:34 PST
- Apparently-to: mathgroup-send at yoda.ncsa.uiuc.edu
>I wonder if it is possible to overlap two surfaces in Mathematica
>(on one plot).
>
>Peter
>flatau at handel.cs.colostate.edu
Basically, the way to do this is to convert both surfaces (SurfaceGraphics
objects) to Graphics3D objects, then Show the resulting Graphics3D objects.
It would be nice to be able to directly combine SurfaceGraphics objects,
but the whole idea of such an object is that if it can be assumed to be a
graph of a SINGLE function of two variables, then the hidden line
algorithms sued in displaying it are much simplified. Try the following:
plot1=Plot3D[Sin[x y],{x,-Pi,Pi},{y,-Pi,Pi}]
plot2=Plot3D[(x^2-y^2)/9,{x,-Pi,Pi},{y,-Pi,Pi}]
plot1=Graphics3D[plot1] (* Converts to Graphics3D object *)
plot2=Graphics3D[plot2]
Show[{plot1,plot2}]
Alternatively, I wrote a widget to do this called MultiPlot3D. It
is in a package called "Student.m", which I will append to this letter.
It is also available for anonymous ftp from otter.stanford.edu in the
/ftp/mma directory.
- Robert Campbell -
Math Dept, UC Berkeley
BeginPackage["Student`"]
(* Copyright 1990 Robert Campbell *)
Bezier::usage = "Bezier[{{x0,y0},...,{x3,y3}}]
computes the cubic bezier curve with control points
{xi,yi}, prints out its parametric description, and
then draws the curve and its control polygon, and
returns the parametric description.
Bezier[{{x0,y0},...},nodraw] returns a list of the two
components of the parametric description."
Midpoint::usage = "Midpoint[f,{x,xmin,xmax},n] will compute the numeric
integral of the function f(x) on the interval [xmin,xmax],
using the midpoint method with n divisions. The numeric
result and the error, or difference from the system command
NIntegrate, will be printed and the result returned will
be a graph of the process. Midpoint[f,{x,xmin,xmax},n,nodraw]
prints the error and returns the numeric result without graphing it."
MultiPlot3D::usage = "MultiPlot3D[{f1,...},{x_,x0_,x1_},{y_,y0_,y1_}]
extends the function Plot3D by allowing the user to plot
several functions of two variables on the same set of axes.
The result is a Graphics3D object, not a SurfaceGraphics object.
The function can use the same options as Plot3D (Not entirely true at
present 14 Feb, 1990). THIS IS A VERY SLOW AND EXPENSIVE PROGRAM."
Newton::usage = "Newton[f,{x,xstart},n] will use n steps
of Newton's method to try to compute a root of f(x),
starting at the point x=xstart and will return a graph of
the process. Newton[f,{x,xstart},n,nodraw] will print the
result at each step of the process, and will return a list
of the values."
Polyapprox::usage = "Polyapprox[f,x,{x0,...,xn}]
will compute the polynomial of degree (n-1) which interpolates
f(x) at the points {x0,...}. It will print this polynomial
and then graph it on the same plot with f along with f."
Trapezoid::usage = "Trapezoid[f,{x,xmin,xmax},n] will
compute the numeric integral of f(x) on the interval
[xmin,xmax], using the trapezoid method with n divisions.
The numeric result and the error, or difference from the
system command NIntegrate, will be printed and a
graph will be returned. Trapezoid[f,{x,xmin,xmax},n,nodraw]
prints the error and returns the numeric result without graphing it."
nodraw;
Begin["`Private`"]
Bezier[{{x0_,y0_},{x1_,y1_},{x2_,y2_},{x3_,y3_}}]:=
(* Compute and display the cubic Bzier curve with the specified control pts.*)
Block[{xcoord,ycoord,plot,points,poly,paramname,t},
(* Define the x and y coordinates of the point on the curve
with parameter un. Generate a new name, of the form un, for the
parametric variable used in Bezier. *)
paramname = Unique["u"];
xcoord=x0*(1-paramname)^3 + 3*x1*paramname*(1-paramname)^2 \
+ 3*x2*paramname^2*(1-paramname) + x3*paramname^3;
ycoord=y0*(1-paramname)^3 + 3*y1*paramname*(1-paramname)^2 \
+ 3*y2*paramname^2*(1-paramname) + y3*paramname^3;
Print["The Bzier curve is:"];
Print["x(",paramname,") = ",xcoord];
Print["y(",paramname,") = ",ycoord];
(* recompute the curve using an internal variable parameter *)
xcoord=x0*(1-t)^3 + 3*x1*t*(1-t)^2 + 3*x2*t^2*(1-t) + x3*t^3;
ycoord=y0*(1-t)^3 + 3*y1*t*(1-t)^2 + 3*y2*t^2*(1-t) + y3*t^3;
(* Now draw the curve invisibly *)
plot=ParametricPlot[{xcoord,ycoord},{t,0,1},DisplayFunction->Identity];
(* Prepare to draw the control points *)
points={Point[{x0,y0}],Point[{x1,y1}],Point[{x2,y2}],Point[{x3,y3}]};
(* Connect the control points to get the control polygon. *)
poly=Line[{{x0,y0},{x1,y1},{x2,y2},{x3,y3}}];
(* combine the graph with the control points, and display them. *)
Show[Graphics[{PointSize[0.02],points,poly}],plot,\
DisplayFunction->$DisplayFunction]
]
Bezier[{{x0_,y0_},{x1_,y1_},{x2_,y2_},{x3_,y3_}},nodraw]:=
(* Compute the cubic Bzier curve with the specified
control pts and return it as a parametric curve in "un".*)
Block[{xcoord,ycoord,paramname},
(* Generate a new name, of the form un, for the parametric
variable used in Bezier. *)
paramname = Unique["u"];
xcoord=x0*(1-paramname)^3 + 3*x1*paramname*(1-paramname)^2 \
+ 3*x2*paramname^2*(1-paramname) + x3*paramname^3;
ycoord=y0*(1-paramname)^3 + 3*y1*paramname*(1-paramname)^2 \
+ 3*y2*paramname^2*(1-paramname) + y3*paramname^3;
{xcoord,ycoord}
]
Midpoint[f_,{x_,xmin_,xmax_},n_]:=
(* This routine will compute the numeric integral of the
function f(x) on the interval [xmin,xmax], using the
midpoint method with n divisions. The numeric result and
the error, or difference from the system NIntegrate, will
be printed and the result returned will be a graph of the
process.*)
Block[{mp,error,i,j,rectlist,linelist,val,step,plot},
mp=0;
step=(xmax-xmin)/n;
For[i=0,i<n,i++,
val[i]=f/. x-> xmin+(i+0.5) step;
mp=mp+val[i] step;
];
rectlist=Table[Rectangle[{xmin+j*step,0},{xmin+(j+1) step,val[j]}]\
,{j,0,n-1}];
linelist=Table[Line[{{xmin+j*step,0},{xmin+j*step,val[j]},\
{xmin+(j+1) step,val[j]},{xmin+(j+1) step,0}}],{j,0,n-1}];
rectlist=Prepend[rectlist,GrayLevel[0.6]];
linelist=Prepend[linelist,GrayLevel[0]];
plot=Plot[f,{x,xmin,xmax},DisplayFunction->Identity];
error=Abs[N[mp-NIntegrate[f,{x,xmin,xmax}]]];
Print["The integral is approximately ",N[mp]];
Print["error=",error];
Show[Graphics[Join[rectlist,linelist]],plot,\
DisplayFunction->$DisplayFunction]
]
Midpoint[f_,{x_,xmin_,xmax_},n_,nodraw]:=
(* This routine will compute the numeric integral of the
function f(x) on the interval [xmin,xmax], using the
midpoint method with n divisions. The error, or
difference from the system NIntegrate, will be printed
and the result returned will be the numeric result. *)
Block[{mp,error,i,j,rectlist,val,step},
mp=0;
step=(xmax-xmin)/n;
For[i=0,i<n,i++,
val[i]=f/. x-> xmin+(i+0.5) step;
mp=mp+val[i] step;
];
error=Abs[N[mp-NIntegrate[f,{x,xmin,xmax}]]];
Print["error=",error];
N[mp]
]
MultiPlot3D[fcns_List,{x_,x0_,x1_},{y_,y0_,y1_},opts___]:=
(* This routine extends the function Plot3D by allowing the
user to plot several functions of two variables on the same
set of axes. The result is a Graphics3D object, not a
SurfaceGraphics object. At present those options which are
valid for Graphics3D, as well as the Mesh option, are not
properly passed or acted on *)
Block[{plotopts,showopts,plotoptstest},
(* Strip out those options which are unique to Plot3D, as
opposed to those options understood by Graphics3D *)
plotoptstest=(MatchQ[#,PlotPoints->_] ||
MatchQ[#,ClipFill->_] || MatchQ[#,HiddenSurface->_] ||
MatchQ[#,Mesh->_] || MatchQ[#,MeshStyle->_])&;
plotopts=Cases[{opts},_?plotoptstest];
showopts=Cases[{opts},_?(Not[plotoptstest[#]]&)];
Show[Map[Graphics3D,Map[Plot3D[#,{x,x0,x1},{y,y0,y1},plotopts,\
DisplayFunction->Identity]&,fcns]],showopts,\
DisplayFunction->$DisplayFunction]
]
Newton[f_,{x_,xstart_},n_,nodraw]:=
(* This routine will use Newton's method to solve the equation f(x)=0 for x.
The initial guess, xstart, is supplied by the user, as is the desired number
of iterations, n. At each step the routine prints out the new estimate, xi, the
value of the function at xi, and the amount the new estimate differs from the
previous estimate. *)
Block[{fprime,xi,xiold,xlist},
xlist = {xstart};
fprime:=D[f,x];
xi=xstart; (* We start out with our first guess *)
For[i=1,i<n,i++,
xiold=xi; (* Keep track of our last guess *)
(* Use Newton's Method to find a better guess.
We keep track of 60 digits internally. *)
xi=xi-N[((f/fprime)/.x->xi),60];
AppendTo[xlist,xi]; (* And add it to our list of guesses *)
Print[i];
Print["xi=",N[xi,10]]; (* Print guess to 10 digits *)
Print["f(xi)=",N[(f /. x->xi),10]];(* Print function value *)
Print["xi-x(i-1)=",N[xi-xiold,10]];(* Change in guess *)
(* Print out the relevant information *)
];
N[xlist,10] (* Return 10 digits as a final value *)
]
Newton[f_,{x_,xstart_},n_]:=
(* This routine will use Newton's method to solve the equation f(x)=0 for x.
The initial guess, xstart, is supplied by the user, as is the desired number
of iterations, n. The output is a graph which shows a plot of the function f,
each successive estimate, xi, and the value of f at xi. The plot also draws
the tangent line to f at {xi,f(xi)} which is used to get the next estimate. *)
Block[{fprime,xi,xiold,i,pointlist,linelist,xmin,xmax,plot},
(* Keep track of max and min x to get the plot bounds *)
xmin = xstart;
xmax = xstart;
fprime := D[f,x];
xi = N[xstart];
(* Put in the initial elements of the graphics lists *)
pointlist={Point[N[{xstart,f/.x->xstart}]]};
linelist={};
For[i=1,i<=n,i++,
xiold = xi; (* Save the old value of xi *)
xi = xi-((f/fprime)/.x->xi);
(* Use Newton's Method to find a better guess *)
(* Draw the line and the new point *)
pointlist=Append[pointlist,Point[{xi,f/.x->xi}]];
linelist=Append[linelist,Line[{{xiold,f/.x->xiold},{xi,0}}]];
(* Keep track of max and min x to get the plot bounds *)
If[xi<xmin,xmin=xi,];
If[xi>xmax,xmax=xi,];
];
(* First we plot the graph invisibly, adding some space on each end *)
plot=Plot[f,{x,xmin-0.2(xmax-xmin),xmax+0.2(xmax-xmin)},\
DisplayFunction->Identity];
(* Then we plot all the desired information *)
Show[Graphics[Join[{PointSize[0.02]},pointlist,linelist]],plot\
,DisplayFunction->$DisplayFunction]
]
Polyapprox[f_,x_,p_]:=
(* This function will find the polynomial of least degree which
passes through the points {x0,f(x0)},{x1,f(x1},... It will print
this polynomial and then plot the polynomial on the same graph
with f(x). Actually, the command Fit can be used to do the same thing.*)
Block[{i,soln,pts,points,a,plot},
pts=Table[{p[[i]],f/.x->p[[i]]},{i,Length[p]}];
poly=Fit[pts,Table[x^i,{i,0,Length[p]-1}],x];
points=Table[Point[{p[[i]],poly/.x->p[[i]]}],{i,1,Length[p]}];
plot=Plot[{f,poly},{x,Min[p]-1,Max[p]+1},DisplayFunction->Identity];
Show[Graphics[{PointSize[0.02],points}],plot,\
DisplayFunction->$DisplayFunction];
poly
]
Trapezoid[f_,{x_,xmin_,xmax_},n_]:=
(* This routine will compute the numeric integral of the
function f(x) on the interval [xmin,xmax], using the
trapezoid method with n divisions. The numeric result and
the error, or difference from the system NIntegrate, will
be printed and the result returned will be a graph of the
process.*)
Block[{trap,error,i,j,rectlist,linelist,val,step,plot},
Array[val,n,0];
mp=0;
step=(xmax-xmin)/n;
val[0]=f/. x->xmin;
For[j=0,j<n,j++,
val[j+1]=f/. x-> xmin+(j+1)*step];
j=.;
trap=Sum[(val[i]+val[i+1])*step/2,{i,0,n-1}];
rectlist=Table[Polygon[{{xmin+j*step,0},{xmin+j*step,val[j]}\
,{xmin+(j+1) step,val[j+1]},{xmin+(j+1) step,0}}],{j,0,n-1}];
linelist=Table[Line[{{xmin+j*step,0},{xmin+j*step,val[j]},\
{xmin+(j+1) step,val[j+1]},{xmin+(j+1) step,0}}],{j,0,n-1}];
rectlist=Prepend[rectlist,GrayLevel[0.6]];
linelist=Prepend[linelist,GrayLevel[0]];
plot=Plot[f,{x,xmin,xmax},DisplayFunction->Identity];
error=Abs[N[trap-NIntegrate[f,{x,xmin,xmax}]]];
Print["The integral is approximately ",N[trap]];
Print["error=",error];
N[trap];
Show[Graphics[Join[rectlist,linelist]],plot,\
DisplayFunction->$DisplayFunction]
]
Trapezoid[f_,{x_,xmin_,xmax_},n_,nodraw]:=
(* This routine will compute the numeric integral of the
function f(x) on the interval [xmin,xmax], using the
trapezoid method with n divisions. The error, or
difference from the system NIntegrate, will be printed
and the result returned will be the numeric result. *)
Block[{trap,error,j,val,step},
Array[val,n,0];
mp=0;
step=(xmax-xmin)/n;
val[0]=f/. x->xmin;
For[j=0,j<n,j++,
val[j+1]=f/. x-> xmin+(j+1)*step];
j=.;
trap=Sum[(val[i]+val[i+1])*step/2,{i,0,n-1}];
error=Abs[N[trap-NIntegrate[f,{x,xmin,xmax}]]];
Print["error=",error];
N[trap]
]
End[]
EndPackage[]