RE: plotting with boundary conditions
- To: mathgroup at smc.vnet.net
- Subject: [mg34511] RE: [mg34482] plotting with boundary conditions
- From: "Wolf, Hartmut" <Hartmut.Wolf at t-systems.com>
- Date: Fri, 24 May 2002 02:42:18 -0400 (EDT)
- Sender: owner-wri-mathgroup at wolfram.com
> -----Original Message----- > From: Tayade Rajeshwary Gunwant [mailto:rajeshwary at neo.tamu.edu] To: mathgroup at smc.vnet.net > Sent: Thursday, May 23, 2002 9:32 AM > Subject: [mg34511] [mg34482] plotting with boundary conditions > > > > Hi all, > I need to plot a surface > f(a,b) = a^2-2*c(a,b)+b^2; given the boundary conditions c(a,b)< ab > > Here c(a,b) is another function of (a,b) given by > c(a,b) = 0.1*(a-1)^2+0.1*(b-1)^2-0.01 > > Is there a direct command to set these boundary > conditions...?? If i give > > f[a_, b_] := If[c[a, b] < a*b, a^2 + b^2 - 2*c[a, b],]; > > Plot3D[f[a, b], {a, 0, 2}, {b, 0, 2}, AxesLabel -> {"a", "b", "f"}, > ViewPoint -> {2, 2, 2}] > > I get error messages for the values of a and b for which > f(a,b) doesnot exist. > > Can anybody please tell me how i can get this plot? > > Thanks > raj > > Raj, I answered to a related question in December last year, see http://library.wolfram.com/mathgroup/archive/2001/Dec/msg00114.html Here, only with few comments, I'll show how that method applies to your case: (1) make a plot disregarding your constraints: In[1]:= f[a_, b_] := a^2 - 2*c[a, b] + b^2 In[2]:= c[a_, b_] := 0.1*(a - 1)^2 + 0.1*(b - 1)^2 - 0.01 In[3]:= g = Plot3D[f[a, b], {a, 0, 2}, {b, 0, 2}] (1a) if you like, look at your constraints: In[4]:= ContourPlot[c[a, b] - a b, {a, 0, 2}, {b, 0, 2}, Contours -> {0}] (2) convert the result to a Graphics3D object: In[5]:= g3D=Graphics3D[g] (3) define three helpers: In[6]:= tagwhere[p : {a_, b_, _}] := If[c[a, b] - a b < 0, {p, inside}, {p, outside}] ...tags points whether inside or outside of your constrained area In[7]:= sol[{p1_, _}, {p2_, _}] := Module[{p, d, a, b}, p = p1 (1 - d) + p2 d; {a, b} = Take[p, 2]; First[ Cases[Solve[c[a, b] - a b == 0, d], s_ /; 0 <= (d /. s) <= 1 :> (p /. s)]]] ...calculates points on the border In[8]:= sep[pp1_, pp2_] := If[pp1[[2]] === pp2[[2]], {pp1}, ppx = sol[pp1, pp2]; {pp1, {ppx, pp1[[2]]}, {ppx, pp2[[2]]}}] ...if two adjacent corner points of a tile are on different sides of the border, calculate the intersection of the connecting edge with the border and add the intersection point to both sides. (4) fix up your graphics and show: In[9]:= gnew1=g3D /. Polygon[pts:{p1_,p2_,p3_,p4_}] :> With[{tagpts = tagwhere /@ pts}, With[{t=sep @@@ Transpose[{tagpts, RotateLeft[tagpts]}]}, With[{newpts= Cases[t, {p_,label_} /; label == inside :> p, 2]}, If[Length[newpts] > 2, Polygon[newpts],{}]]]] In[10]:= Show[gnew1] -- Hartmut