title: How can I effectively performance-profile Dynamics?

text: How to optimize performance-wise a complex dynamic expression (like a GUI) with lots of dynamically updated components, which require different update- frequencies? I was hesitant to ask this question as it does not focus on a very specific problem, but is about general speed (and memory) improvement tricks. The gui is a primitive scrolling racer where the user has to manipulate a controller to keep an object on a winding path. I did it for a cognitive experiment for humans but it definitely has potential in the gaming market :) Also - as far as I know - _Mathematica_ was not used (and was not designed) to develop realtime platformer/racer games (or alike), so the question naturally arises: can _Mathematica_ do it? This, I think, is untrodden land worthy of exploration. Below is a toy model of a game. I tried to keep it as abstract as possible to make it general and thus applicable to a wider domain. If you still think that this question is too localized/not really a question, please say so and I remove it. **Details of the code** The code might seem too long, but I assure you, it is not complicated. If the game is started with the Start button, the function `start` initiates a `ScheduledTask` that sets the in-game timer with `updateFrequency = .05` seconds. The controller can be dragged continuously, and it generates a coordinate pair ${x_1, x_2}$, stored in `ctrl`, that is fed to the `controlFunction`, that defines the actual vertical position (`y`) of the yellow arrow (the "ship"). The responsible equation is: $$ y = y* + k_1 x_1 + k_2 x_2 $$ where $y*$ is the $y$ value of the position of the arrow in the previous timestep. Now the hard part is, that from time to time, parameters $k_1$ and $k_2$ are randomly changed (every 100 steps, using the `change` function), so the controller is changed as well. The controller background always displayes the actual gradient of ${\varDelta}y$. The different background colors in each episode indicate that the controller has changed. I did not use `Locator` or `LocatorPane`, as those seemed to be unresponsive, especially when the underlying controller is changed. While the controller is updated continuously, the actual position of the arrow is only updated via the `update` function when the running `ScheduledTask` tells so (in every 0.05 sec). Also, moving the controller produces some sounds that helps the user to learn the control gradient.               DynamicModule[{       optimumFunction = Sin[#*\[Pi]/100] &, (* defines the winding optimum zone *)       updateFrequency = .05, (* update plot every 0.05 seconds *)       xFrame = 100, (* horizontal range to show (last 100 steps) *)       yFrame = 30, (* default vertical range to show *)       maxTime = 1000, (* length of a run. Real word time is ~ 0.05 * 1000 = 50 sec *)       episodeNum = 5, (* number of different episodes: each episode lasts 200 timesteps *)            active, controlActive, ctrl, data, time, \[Theta], \[Rho], y,        yRange, xRange, gradient, crosshair, ship, episodeColors,        episodeBorders,       reset, sound, start, stop, update, change, controlFunction       },           (* Internal functions *)      reset[] := (active = False; ctrl = {0, 0}; data = {}; time = 0; y = 0; change[]; update@0);      sound[{x_, y_}, vol_] := Sound[SoundNote[Floor@Rescale[controlFunction[x, y], {-1, 1}, {-20, 40}], .3, "Strings"], SoundVolume -> vol];      start[] := (active = True; RunScheduledTask[update@(time++), {updateFrequency, \[Infinity]}]);      stop[] := (RemoveScheduledTask[ScheduledTasks[]]; active = False);      update[t_] := Block[{xAll, yAll, n = 2},        y = y + controlFunction @@ ctrl; (* update the position of the ship *)        data = Append[data, {t, y}]; (* save position of the ship *)        yRange = ({Min[Min@#, -{-1.5, 1.5}], Max[Max@#, {-1.5, 1.5}]}*1.3) &@(Last /@ Take[data, -Min[yFrame, Length@data]]); (* calculate vertical and horizontal ranges to show on the scrolling plot*)        xRange = {Max[0, time - xFrame], Max[xFrame, time]};        xAll = N@Abs@(Subtract @@ xRange); (* absolute x-range shown *)        yAll = N@Abs@(Subtract @@ yRange); (* absolute y-range shown *)        \[Theta] = Mod[ArcTan @@ #, 2 \[Pi], 0] &@ If[t < n, {1, 0}, (data[[-1]] - data[[-n]])/{xAll, yAll}]; (* rotation angle of the arrow, dependent on previous steps and the actual plot range used *)        If[0 < t < maxTime \[And] Mod[t, maxTime/episodeNum] == 0, change[]]; (* change controllers *)        If[t >= maxTime, stop[]]; (* stop run *)        ];      change[] := (        \[Rho] = RandomReal[{-.5, .5}, 2]; (* \[Rho] = {k1, k2}, the two random parameters of the gradient *)        controlFunction = {#1, #2}.\[Rho] &; (* the control function defining the gradient *)        gradient =          DensityPlot[(controlFunction[x, y]), {x, -1, 1}, {y, -1, 1},           Frame -> False, PlotRangePadding -> 0, ImagePadding -> 0,          ImageSize -> 200,           ColorFunction -> ((Lighter@               ColorData["Rainbow"][Rescale[#, {-1, 1}]]) &),           ColorFunctionScaling -> False,           PlotRange -> {{-1, 1}, {-1, 1}, {-1, 1}}, Mesh -> {{{0, Thick}}},           MeshFunctions -> {#3 &}];        );                (* Initialization *)      crosshair =        Graphics[{White, Disk[{0, 0}, 1.5], Darker@Blue, Disk[], Blue,          AbsoluteThickness@1, Circle[{0, 0}, 1.5], Line@{{0, -2}, {0, 2}},          Line@{{-2, 0}, {2, 0}}}, ImageSize -> 15];      ship = Graphics[{FaceForm@Yellow, EdgeForm@{Gray, Thick}, Polygon@{{0, -1}, {4, 0}, {0, 1}, {0, -1}}}, ImageSize -> 30];      episodeBorders = Partition[Range[0, maxTime, maxTime/episodeNum], 2, 1];      episodeColors = {Hue[.66, 1, 1], Hue[.13, .8, 1], Hue[.23, 1, .8], Hue[.5, .4, .6], Hue[.74, 1, .8]};      reset[]; (* initialize first episode *)                (* Output *)      Panel@Grid[{{          Column@{                 (* Buttons *)            Row@{              Dynamic@Button[If[active, "Stop", "Start"], If[active, stop[], start[]]],              Button["Reset", stop[]; reset[]]              },                 (* Controls *)            Framed@Dynamic[              If[! active,                Graphics[{}, Frame -> False, ImageSize -> 200,                 ImagePadding -> 5, Background -> GrayLevel@.9],               EventHandler[Dynamic[                 If[controlActive, ctrl = Clip@MousePosition@"Graphics"; EmitSound@sound[ctrl, .5]];                 Graphics[{                   Texture@gradient,                    Polygon[{{-1, -1}, {-1, 1}, {1, 1}, {1, -1}},                     VertexTextureCoordinates -> {{0, 0}, {0, 1}, {1, 1}, {1, 0}}],                   Dynamic@Inset[crosshair, ctrl]                   }, Frame -> False, ImageSize -> 200, ImagePadding -> 0],                 TrackedSymbols :> {controlActive}],                {                 "MouseDown" :> (controlActive = True;),                  "MouseUp" :> (controlActive = False;)                }]],              TrackedSymbols :> {active, gradient}]            },               (* Scrolling pane *)          Panel[Dynamic[            Show[             (* background colors *)             Graphics[MapThread[               {FaceForm@{Opacity@.2, #1}, EdgeForm@None,                  Rectangle @@ {{First@#2, -100}, {Last@#2, 100}}} &,               {episodeColors, episodeBorders}],               PlotRange -> {Automatic, yRange}],                  (* initial/terminal gray edge (before time 0 and after maxTime *)             Plot[0, {t, 0 - maxTime*.1, 0}, PlotStyle -> None, FillingStyle -> GrayLevel@.9, Filling -> {1 -> Top, 1 -> Bottom}, PlotRange -> {Automatic, yRange}],             Plot[0, {t, maxTime, maxTime*1.1}, PlotStyle -> None, FillingStyle -> GrayLevel@.9, Filling -> {1 -> Top, 1 -> Bottom}, PlotRange -> {Automatic, yRange}],                  (* desired zone (winding path) *)             Plot[{optimumFunction@t - .5, optimumFunction@t,                optimumFunction@t + .5}, {t, Max[0, time - xFrame] - 20,                Max[xFrame, time] + 100},               PlotStyle -> {None,                 Directive[GrayLevel[0, .7], AbsoluteThickness@3, Dashed],                 None}, FillingStyle -> Hue[.3, 8, .8, .5],               Filling -> {1 -> {3}}],                  (* ship & trajectory *)             Graphics[{GrayLevel@.6, AbsoluteThickness@2,                Dynamic@Line@Most@data,                Inset[Rotate[ship, \[Theta]], Last@data]}],                  (* options *)             PlotRange -> {xRange, yRange},              PlotRangePadding -> {{Scaled@.02, Scaled@.2}, 0},             ImageSize -> Full, AspectRatio -> Full,              ImagePadding -> Automatic,             Frame -> True, FrameTicks -> False, Axes -> False,              GridLines -> {Range[0, maxTime, maxTime/episodeNum], None},              Method -> {"GridLinesInFront" -> True}             ],             TrackedSymbols :> {yRange, xRange, time, y, \[Theta], active}            ], ImageSize -> {800, 400}]          }}, Alignment -> {Left, Center}]]      ![Mathematica graphics](http://i.stack.imgur.com/DmC4C.png) I know that the best advice here is that   * only update if necessary;   * only update the smallest minimum that needs to be updated. I was keeping that in mind during development. Most of the scrolling graphics are pre-generated (the sine wave and background colors as rectangles) and only the `PlotRange` is changed dynamically. Explicit `TrackedSymbols` are used to update only those symbols that should be updated. Still, the gui is lagging: the controller follows the moused with a reasonably large lag, that does not feel smooth. Also, during every episode change, there is a short pause. Now, this lagging is inacceptable: my first Commodore had only 64 KB memory, but it could run such games way smoother than _Mathematica_.   1. **How can I understand what uses up time during dynamic evaluation (considering that I do not want to print`AbsoluteTiming` values all the time)? How can I profile such a gui?**   2. **How can I improve the code further to give a smooth feeling?** Any comments/ideas are welcome. * * * I should tell here that this code uses a lot of stuff from all around Mathematica.SE, so if you run into some familiar lines, don't be surprised. Thanks for all your helpful input. Also thanks for Oleksandr for the title! * * * **EDIT at 7. 8. 2012** According to Rojo's profiling and comments in chat, here is a modified version of the scrolling output (just replace the original internal `Panel[Dynamic[Show[...]]]` with the code below. It removes the wrapping `Dynamic` around `Show` and instead uses a dynamically updated `PlotRange` (also I moved the `Dynamic`s in the ship and trajectory to the inside of `Inset` and `Line`). Now while this updates much faster, and any check on `xRange` and `yRange` indicates that they are updated correctly at every timestep, the displayed plot range is only updated when the controller is moved, meaning that if the controller is not touched, the ship will just move out of the range on the right (clearly unwanted).               Panel[Show[       Graphics[MapThread[         {FaceForm@{Opacity@.2, #1}, EdgeForm@None, Rectangle @@ {{First@#2, -100}, {Last@#2, 100}}} &,         {episodeColors, episodeBorders}]],            (* initial/terminal gray edge *)       Plot[0, {t, 0 - maxTime*.1, 0}, PlotStyle -> None, FillingStyle -> GrayLevel@.9, Filling -> {1 -> 100, 1 -> -100}],       Plot[0, {t, maxTime, maxTime*1.1}, PlotStyle -> None, FillingStyle -> GrayLevel@.9, Filling -> {1 -> 100, 1 -> -100}],            (* desired zone ("road") *)       Plot[{optimumFunction@t - .5, optimumFunction@t, optimumFunction@t + .5}, {t, 0, maxTime}, PlotStyle -> {None, Directive[GrayLevel[0, .7], AbsoluteThickness@3, Dashed], None}, FillingStyle -> Hue[.3, 8, .8, .5], Filling -> {1 -> {3}}],            (* ship & trajectory *)       Graphics[{GrayLevel@.6, AbsoluteThickness@2, Line@Dynamic@Most@data,Inset[Rotate[ship, Dynamic@\[Theta]], Dynamic@Last@data]}],            (* options *)       PlotRange -> Dynamic[{xRange, yRange}, TrackedSymbols :> {xRange, yRange}],        PlotRangePadding -> {{Scaled@.02, Scaled@.2}, 0},       ImageSize -> Full, AspectRatio -> Full, ImagePadding -> Automatic,       Frame -> True, FrameTicks -> False, Axes -> False,        GridLines -> {Range[0, maxTime, maxTime/episodeNum], None},        Method -> {"GridLinesInFront" -> True}       ], ImageSize -> {800, 400}]


title: How to preserve the order of terms in the output

text: Is there a way in which mathematica can preserve the order of operations while it evaluates an expression?, for instance: In my short example               w = f[r,Θ] ;        lapla12 = -(r/(2 μ)) (1/r D[w/r r, {r, 2}] - 1/r^2 (D[w/r, {Θ, 2}] + Cot[Θ] D[w /r, {Θ, 1}])) // ExpandAll // TraditionalForm      The output that mathematica gave me is ![enter image description here](http://i.stack.imgur.com/M6h9i.jpg) However, I know the real order is ![enter image description here](http://i.stack.imgur.com/kdznj.jpg) In my real problem, I have more than 100 terms and is quite difficult to me sort all terms to check the results. **EDITED:** For instance the first 35 terms in my real expression appears in this order and i would like that they appears in the correct order![enter image description here](http://i.stack.imgur.com/V9EeQ.jpg) I'm using the code from Jens for derivatives formatting purposes. Here is the code i'm using to obtain the above output.               ClearAll["Global`*"]     SetOptions[$FrontEndSession, PrintingStyleEnvironment -> "Condensed"]         $PrePrint = # /. {Csc[z_] :> 1/Defer@Sin[z],           Sec[z_] :> 1/Defer@Cos[z],           Cot[z_] :> Defer@Cos[z]/Defer@Sin[z]} &;               Subscript[w, 1] = f[Subscript[r, 1], Subscript[r, 2], α, β, γ, Subscript[Θ, 12]];          Set[J, -(D[Subscript[w, 1], {β, 2}] + 1/Sin[β]^2 (D[Subscript[w, 1], {α, 2}] +               D[Subscript[w, 1], {γ, 2}]) + Cot[β] D[Subscript[w, 1], β] -            2 Cot[β]/ Sin[β] D[Subscript[w, 1], α, γ])] // TraditionalForm ;          Set[Subscript[a, 1], 2 D[Subscript[w, 1], {β, 2}] + D[Subscript[w, 1], {γ, 2}] + J]          // TraditionalForm ;          Set[Subscript[a, 2],  2 Cot[β]/ Sin[β] D[Subscript[w, 1], {α, 1}] - (1 + 2 Cot[β]^2)          D[Subscript[w, 1], {γ, 1}] - 2/Sin[β] D[Subscript[w, 1], β, α] + 2 Cot[β]          D[Subscript[w, 1], β, γ]] // TraditionalForm ;          Set[Subscript[f, 1], 1/(2 Sin[Subscript[Θ, 12]]^2) (-J + Cos[Subscript[Θ, 12]]          (Sin[2 γ] Subscript[a, 2] - Cos[2 γ] Subscript[a, 1]) + Sin[Subscript[Θ, 12]]         (Sin[2 γ] Subscript[a, 1] + Cos[2 γ] Subscript[a, 2])) - D[Subscript[w, 1],          Subscript[Θ, 12], γ] - 1/2 Cot[Subscript[Θ, 12]] D[Subscript[w, 1], {γ, 1}] +         (1/4 - 1/(2 Sin[Subscript[Θ, 12]]^2)) D[Subscript[w, 1], {γ, 2}]          ] // ExpandAll // TraditionalForm ;          Set[lapla1, -(1/(2 Subscript[μ, 1])) (1/Subscript[r, 1] D[Subscript[w, 1]          Subscript[r, 1], {Subscript[r, 1], 2}] - 1/Subscript[r, 1]^2         (D[Subscript[w,  1], {Subscript[Θ, 12], 2}] + Cot[Subscript[Θ, 12]]          D[Subscript[w, 1], {Subscript[Θ, 12], 1}]) - Subscript[f, 1]/Subscript[r, 1]^2)         ] // ExpandAll // TraditionalForm      **EDITED 2** In this case the correct order is ruled by the expression for lapla1.![enter image description here](http://i.stack.imgur.com/TRGS1.png) where terms in red (radial terms) go first, then terms in purple (internal angular terms), and finally the terms in blue (coupling between the orientational angular terms and the internar angular terms). Blue terms are ruled by the expression for $f_1$ from left to right.


title: NDSolve: methods and step size choosing

text: I am looking into the documentation of `NDSolve[]`; more precisely **how this function chooses the`StepSize` and how it chooses which method to use**. * * * # Regarding `StepSize` control, it reads > NDSolve adapts its step size so that the estimated error in the solution is > just within the tolerances specified by PrecisionGoal and AccuracyGoal. **Question:** How? Does it use an adaptive step size similar to this one for instance? * * * # Regarding method choice, it reads > `NDSolve` typically solves differential equations by going through several > different stages depending on the type of equations. With `Method->{s1->m1, > s2->m2, ...}` stage is handled by method . The actual stages used and their > order are determined by `NDSolve`, based on the problem to solve. and also > Possible explicit time integration settings for the `Method` option include: > [Adams, BDF, etc] with various submethods and options. A good tutorial about method choice and stiffness detection is found here. **Questions:**   1. What is the default method choice in `NDSolve`? (I am guessing it will start by `ExplicitRungeKutta`, seeing as it is the most robust from what I've been reading, but then which?)   2. Regarding stiffness detection and method choice upon stiffness finding, I found this page that refers "Explicit modified midpoint, double-harmonic sequence, linearly implicit Euler, and sub-harmonic sequence" as possible methods. In that same page, I found how those methods seem to compute stability of the model, and use it to determine whether it should use a stiff/non-stiff method. Yet, what is again the default order of choice for those methods? Thank you for your insight.


title: J/Link: How to run a thread inside Mathematica

text: In the J/Link User Guide there is no example on how to use java.lang.Thread, from inside Mathematica. I'd like to know if is possible to start a thread from a `JavaNew` object which has a Runnable interface implemented. If I run this code from a java main class everything works fine. The Thread will run sampling bid/ask from a tcp/ip data stream and every time I need to know the bid/ask data I call the method that get me the value.               t3OpenT job01 = new t3OpenT("192.168.0.78",5333,500000,200);             job01.openConnection();             job01.setRequest("function=subscribe|item=MI.EQCON.1|schema=best_bsiz1;best_bid1;best_ask1;best_asiz1;symbol");             job01.subscribeData();          System.out.println("now start the thread 01");             Thread ntJob01 = new Thread(job01);             ntJob01.setDaemon(true);             ntJob01.setName("GENERALI");             ntJob01.start();      But if I create the objects inside the notebook using JavaNew[....] and then start the thread it will execute but it is interrupted immediately. Here is the mathematica code I use:               In[33]:= generaliGRT = t3OpenObjectT[{generaliG}]          Out[33]= {{"1", "G", JLink`Objects`vm1`JavaObject27681839799861249}}          (*Now we setup T3OpenObject to subscription:*)          In[34]:= generaliGRT[[1, 3]]@      setRequest[       "function=subscribe|item=MI.EQCON.1|schema=best_bsiz1;best_bid1;best_ask1;\     best_asiz1;symbol"]          (*subscribe the data on GENERALI stock*)          In[35]:= generaliGRT[[1, 3]]@subscribeData[];     generaliGRT[[1, 3]]@getResponse[]          Out[36]= "outcome=OK|item=MI.EQCON.1"          (*the response is OK so now we are redy to start the thread *)          In[37]:= th = JavaNew["java.lang.Thread", generaliGRT[[1, 3]]]          Out[37]= JLink`Objects`vm1`JavaObject8825827030466561          th@setDaemon[True]          In[39]:= th@setName["GENERALI"]          In[40]:= th@getName[]          Out[40]= "GENERALI"          In[41]:= th@start[]          In[42]:= th@getId[]          Out[42]= 14          (*the thread is immediately terminated *)     In[43]:= th@getState[]@toString[]          Out[43]= "TERMINATED"          (*but I can read the pusched stock date*)      In[45]:= generaliGRT[[1, 3]]@getInfo[]          Out[45]= {"MI.EQCON.1", "18387", "17.19", "17.21", "14300", "G"}      and here is the run() method of T3OpenT class:               @Override         public void run(){                  thExtSignal = false;             refThread = Thread.currentThread();                  try{                      while (true){                          try{                         //check if thread must be terminated                         if(thExtSignal == true){                             Thread.currentThread().interrupt();                         }                         if(Thread.interrupted()){                              throw new InterruptedException();                         }                         //read data from buffer                         this.pushedData = this.refBufReader.readLine();                         //get the time stamp                         dateTimeSample = dateFormat.format(cal.getTime()).toString();                         //split the response end move data in a string vector                         String responseBidAsk[] = pushedData.split("\\|", 6);                         //update only elements tha are changed                         for(int i=0 ; i<=responseBidAsk.length-1;i++){                             if(responseBidAsk[i].length()>0){                                 info[i]=responseBidAsk[i];                             }else{                                 info[i]=info[i];                             }                         }                     }catch (InterruptedException e) {                         //if the thread is interrupted then unsubscribe                         this.unsubscribe();                         // We've been interrupted: no more messages.                         return;                             }                 }             }                catch(SocketTimeoutException e){                 //if timeout unsubscribe                 response = "timeOut";                 unsubscribe();             }             catch (IOException e) {                 pushedData = "error";                  }              }


title: How to do DateDifference on table items containing dates

text: I have a `Table` with data in form `{{date}, value}`:               data= {{{2012, 6, 1}, 16}, {{2012, 6, 8}, 14.24}, {{2012, 6, 15}, 13.7},             {{2012, 6, 22}, 12.31}, {{2012, 6, 29}, 11.5}, {{2012, 7, 6}, 10.08},             {{2012, 7, 13}, 9.18}, {{2012, 7, 20}, 8.65}, {{2012, 7, 27}, 7.8},             {{2012, 8, 3}, 7.51}, {{2012, 8, 10}, 7.31}, {{2012, 8, 17}, 6.92},             {{2012, 8, 24}, 6.57}, {{2012, 8, 31}, 6.1}, {{2012, 9, 7}, 5.67},             {{2012, 9, 14}, 5.54}, {{2012, 9, 21}, 5.17}, {{2012, 9, 28}, 4.83},             {{2012, 10, 5}, 4.62}, {{2012, 10, 9}, 4.53}}      How can I know the difference (in days) between each date and the starting date `{2012, 6, 1}` and associate this new value to each row of the `Table`?


title: Plot 2D Vector function in 3D

text: I want to plot a 2D vector function such as $F(x,y) = (a(x,y),\,b(x,y))$ in a 3D graph so that the vectors are embedded in the xy plane. I tried to do the following: First I defined a piecewise function like this               g[z_] := Piecewise[{{1, z == 0}}, 0]      Then I converted the 2D vector function to a 3D one by setting the 3rd component to zero and multiplying 1st and 2nd components by `g[z]` so that the x and y components are null when `z != 0`:               f[x_, y_, z_] := {x g[z], y g[z], 0};      Plot the function:               VectorPlot3D[{x g[z], y g[z],0}, {x, -5, 5}, {y, -5, 5}, {z, -5, 5}]      The issue with this solution is that `VectorPlot3D` won't evaluate the function in the relevant points, The above example shows an empty graph, because _Mathematica_ jumps from `z = -1` to `z = 1` without evaluating `z = 0`. I tried with `RegionFunction` (which would've rendered the definition of the above-mentioned piecewise function useless), but that only accepts inequalities, and I want to evaluate the function at any coordinate `{x ,y, 0}`. I could feed it a list of vectors via `VectorPoints -> { {a, b, 0}, {c, d, 0}, ...}`, but that's not an elegant solution at all. Are there other ways to do this?


title: Unwanted PlotMarkers appear at non-existent places

text: In this simple code               ListPlot[Table[Sin[x], {x, 0, Pi, .1}], PlotRange -> {0, 0.5},      Joined -> True, PlotMarkers -> Automatic, ClippingStyle -> None, Frame -> True]      ![enter image description here](http://i.stack.imgur.com/r2nbW.png) We can see that there are two PlotMarkers appear at the clipped boundary. But actually there are not points, so they should not exist. How can we remove them? I have already used `ClippingStyle -> None`. * * * **Update** As Alexey Popkov pointed out, this is really a bug. If we play with this code               data = {Table[{x, Sin[x]}, {x, 0, 1.05 Pi, .1}], Table[{x, Sin[2 x]}, {x, 0, 1.05 Pi, .1}]};     ListPlot[data, PlotRange -> {0, 0.5}, Joined -> True, PlotMarkers -> {Automatic, Medium}]      ![enter image description here](http://i.stack.imgur.com/G7beA.png) Now, even the same plot marker will be used for all data sets at the clipped boundaries, both top and bottom.


title: How can I get thicker lines in the legends of a ListLinePlot?

text: I want to have thicker lines in the legends of this plot:               historyobjetivomejor = {10, 15, 20, 30, 35}     historyobjetivomedia = {6, 8, 10, 13, 15}     historyobjetivopeor = {6, 3, 5, 7, 8}          ListLinePlot[{historyobjetivomejor, historyobjetivomedia, historyobjetivopeor},        PlotStyle -> {Thickness[0.008]},        AxesLabel -> {"Generacion", "Objetivo"},        LabelStyle -> {24, Bold},        ImageSize -> Large,        PlotLegends -> {"Max", "Mea", "Min"},        PlotStyle -> Thickness[0.005],        AxesStyle -> Thickness[0.005]]      I was trying with options from other posts, but without good result My last attemp was               ListLinePlot[{historyobjetivomejor, historyobjetivomedia, historyobjetivopeor},         PlotStyle -> {Thickness[0.008]},        AxesLabel -> {"Generacion", "Objetivo"},        ImageSize -> Large,        PlotLegends -> [{"Max", "Mea", "Min"},        LabelStyle -> {Bold, 20},        LegendMarkerSize -> {40, 10}],        PlotStyle -> Thickness[0.005],        AxesStyle -> Thickness[0.005]]


title: Assign a Distribution to a group of symbolic variable

text: I'm looking to compute this Expectation : ![enter image description here](http://i.stack.imgur.com/5xKXu.png) "i" represent a individual and h(i) is value. I would like the compute the Expected value of his value times the mean of the groups. Everything is unknown so symbolic is conserved. Any help would be greatly appreciated Thank Alexandre


title: Clustering of peaks in 2D spectra

text: I'm having some problems clustering data in a 2D spectra and will really appreciate some help. I've been playing around with the FindClusters function but wasn't able to tune it appropriately. Maybe someone can give me a hand on how to customize the DistanceFunction option or so. A sample of my data looks like this: ![Sample data](http://i.stack.imgur.com/7J8ov.jpg) Ideally I would want the peaks to be clustered as this: ![Ideal clustering](http://i.stack.imgur.com/Y5JFW.jpg) I tried the FindClusters functions with different parameters but only got partial clusters like these: ![Partial clusters](http://i.stack.imgur.com/Sm3mJ.jpg) I believe it may be easy if I can define my own DistanceFunction but I don't understand how to use a custom function there with a threshold or so. I've seen other posts like this one, this other one and this one too but I wasn't able to successfully apply those methods. It doesn't matter if I don't get that red circle clustering like the ideal image. If I can just get to identify the whole clusters, I can reduce that set to just one point that would be the center of gravity. It would be awesome if you can point me the way on how to solve this. You can download that sample from here. Just one thing, I don't know before hand how many clusters I'll get. I believe there are 30 in this dataset but this is just a sample. Thank you!!


title: Removing Abs from Abs[a + Exp[I*c]b]^2

text: Title says it all. Mathematica produces sometimes expressions of the form               Abs[a + Exp[I*c]b]^2      where all three quantities a, b and c are positive real numbers. That expression of course reduces simply to the law of cosines               a^2 + b^2 + 2a b Cos[c]      however, I cannot get Mathematica to display it in this form. Any ideas?


title: How do I force all Graphics Text to be StandardForm?

text: I've been playing around for while trying to determine how to get all graphic objects to use StandardForm fonts. That is, I want all axes labelling, titling, legends, inserted text, etc., for 2D/3D graphics to use StandardForm. Forever and always unless explicitly overridden. I'm getting somewhat overwhelmed with all of the options in the Option Inspector and if anyone has a suggestion I would appreciate it.


title: Trouble with Minimize

text: Minimize[{t^2, 1 == 3*t + 16*u}, {t, u}, Integers]         1. It does not give the solution t = -5 and u = 1.   2. Am I missing something? The solution being known and easy, the actual question is not "how can I minimize?"; it is rather "how is it that _Mathematica_ , even version 8 or 9, can't do it?"


title: Adjust the results of evaluating Table so it can be used as a argument to Point

text: I have a list generated by `Table` as below.               Table[{a x1 + b x2 + c x3, a y1 + b y2 + c y3, a z1 + b z2 + c z3},        {a, -5, 5}, {b, -5, 5}, {c, -5, 5}]      but it can't be used as an argument to `Point`. That is,               Point[       Table[{a x1 + b x2 + c x3, a y1 + b y2 + c y3, a z1 + b z2 + c z3},          {a, -5, 5}, {b, -5, 5}, {c, -5, 5}]]      doesn't work. How do I adjust the list returned by `Table` to make compatible with `Point`?


title: How can I compute and plot the 95% confidence bands for a fitted logistic regression model?

text: I have a set of binary data that I' d like to fit a logistic regression model to and show the 95% confidence band around the fitted function. I can produce the regression model but how do I compute and show the confidence bands in plot3 ? If showing the confidence bands is not possible can the single prediction bands be shown ? Many thanks.               damagecodedata = {{53, 1}, {57, 1}, {58, 1}, {63, 1}, {66, 0}, {67, 0},                       {67, 0}, {67, 0}, {68, 0}, {69, 0}, {70, 1}, {70, 0},                       {70, 1}, {70, 0}, {72, 0}, {73, 0}, {75, 1}, {75, 0},                       {76, 0}, {76, 0}, {78, 0}, {79, 0}, {81, 0}};          plot1 = ListPlot[damagecodedata, PlotRange -> {{50, 90}, {-.1, 1.1}},                Frame -> True,               FrameStyle -> Directive[FontSize -> 13],                FrameLabel -> {Style["Temperature (deg F)", 14], Style["Coded Damage", 14]},                FrameTicks -> {Range[50, 90, 5], {0, 1}, None, None},               PlotStyle -> {PointSize[Medium], Red}, ImageSize -> 400             ]          model = LogitModelFit[damagecodedata, x, x]     Normal[model]     model["ParameterTable"]          plot2 = Plot[model[x], {x, 50, 90}];     plot3 = Show[plot1, plot2]      ![enter image description here](http://i.stack.imgur.com/rk82k.jpg)


title: Poor Graphics3D rendering

text: I am rendering a world environment simulation in _Mathematica_ with `Graphics3D` code:               Show[{     Graphics3D[{     Lighting -> "Neutral",     Green,      Polygon[{{-10, -10, 0}, {-10, 10, 0}, {10, 10, 0}, {10, -10,         0}}],     Red, Cone[{{0, 0, 0}, {0, 0, .1}}, 1/2],     White,      Table[Sphere[        Table[{0, 0, 10} + {rx, ry, rz} + {r, r, r/3} /.           r :> RandomReal[]*4, {5}],         2] /. {rx -> RandomReal[{-20, 20}],         ry -> RandomReal[{-20, 20}], rz -> RandomReal[{-1, 3}]}, {30}]     }, Boxed -> False, Background -> Blue,      ImageSize -> {900, 900*41/54}, SphericalRegion -> True],     ContourPlot3D[     x^2 + y^2 == 10^2, {x, -10, 10}, {y, -10, 10}, {z, 0, .6},      Mesh -> None, Lighting -> "Neutral",      ContourStyle -> Darker[Green]]     }]      If the camera is positioned far away from the scene so that everything is in view range, it renders fine: ![somewhat cloudy](http://i.stack.imgur.com/rtbM4.png) However, when I specify the desired camera parameters:               ViewVector -> {{0, 0, .1}, {0, 1, .1}}, ViewAngle -> 54 \[Degree]      The rendering becomes weird: ![weird rendering](http://i.stack.imgur.com/q02aw.png) The bottom part of the green rectangle is cut off and the blue background is displayed instead, and the red cone is very weird as well. Is there a way to display this scene properly? edit: **Solved** with               ViewRange -> {0.1, 100}      This is not optimal, because it does not take `0` or `Automatic` for a minimum value, `All` causes the cut-off in the above example, and it does not take `Infinity` or `Automatic` for a maximum value.


title: Significant platform-dependent efficiency when importing data?

text: I'm curious, as I understand that answer this way: > If you are using a Linux derivative or something like cygwin on Windows, > then this is an efficient solution for reading only the third field from > your files ( assumes whitespace separated data, other delimiters can be > specified ).               data= ReadList["!head --lines=-2 " <>#<> " | tail --lines=+39 | cut -f 3 "]& /@ FileNames[]      If one deals with data files in the 2 digit GB range, and often only some columns of the data are needed for further evaluation in experimental science, should one think about running mathematica on Linux, cygwin,... Is the efficiency boost for importing data small, big, huge? Do I have more options to select specific data elements in files?


title: Solving a system of equations

text: Mathematica doesn't want to solve my exact system with 3 equations and 3 variables.               AdjustPar20[x_] := a x^b Exp[c x]     Solve[AdjustPar20[0.05] == 9.37126 && AdjustPar20'[0.05] == 0 &&        AdjustPar20[0] == 0, {a, b, c}]      As an explanations why this cannot be solved the program replies:               {}     Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found;                   use Reduce for complete solution information. >>      Any ideas?


title: How to collect terms with positive powers in polynomial

text: I am trying to collect all terms with non-negative powers of $x$ in polynomials like $\frac{1}{x^2}\left(a x^2+x^{\pi }+x+z\right)^2$ First expand the polynomial               Expand[1/x^2 (x + x^π + a x^2 + z)^2, x]      This gives $a^2 x^2+2 a x^{\pi }+2 a x+2 a z+2 x^{\pi -2} z+x^{2 \pi -2}+2 x^{\pi -1}+\frac{z^2}{x^2}+\frac{2 z}{x}+1$. Now try to select those terms with positive or zero powers of x. My best guess is               Plus@@Cases[%, (Times[___, Power[x, g_.], ___] /; g >= 0)]      However, this only yields the term $a^2 x^2$. And why does this not work for the other terms? How can I collect the other terms with positive or zero powers of x?


title: how to superimpose two best-fit graphs

text: I have a couple of best fit curves. All I want to do is superimpose them, to see if they are shaped differently. (I can eyeball them but I want to be sure, as it were). On to what I am asking: I have a couple of data sets. I graphed them both with the following and all was well:               IntData = {{0.05`, 2.835`}, {0.1`, 2.84`}, {0.2`, 2.845`}, {0.3`, 2.85`},                {0.4`, 2.865`}, {0.5`, 2.88`}, {0.6`, 2.9`}};     Fit[IntData, {1, amp, amp^2}, amp]     Plot[2.84 - 0.00005*amp + 0.176*amp^2, {amp, 0, 0.6},       AxesLabel -> {"Amplitude", "Period"}, AxesOrigin -> {0, 2.84}]      I did it again with a second data set. Everything looks the same except the numbers are different: ˚ IntData = Table[{amplitude[n], period[n]}, {n, 0, 6}]; Fit[IntData, {1, amp, amp^2}, amp] Plot[4.02 - 0.00698*amp + 0.271*amp^2, {amp, 0, 0.6}, AxesLabel -> {"Amplitude", "Period"}, AxesOrigin -> {0, 4.02}] Both sets generated nice curves. I wanted to have the curves superimposed so I tried setting up both data sets again with variables named differently:               amplitude1[0] = 0.05; period1[0] = 2.835;     amplitude1[1] = 0.10; period1[1] = 2.840;     amplitude1[2] = 0.20; period1[2] = 2.845;     amplitude1[3] = 0.30; period1[3] = 2.850;     amplitude1[4] = 0.40; period1[4] = 2.865;     amplitude1[5] = 0.50; period1[5] = 2.880;     amplitude1[6] = 0.60; period1[6] = 2.900;     IntData = Table[{amplitude1[n], period1[n]}, {n, 0, 6}];     Fit[IntData, {1, amp1, amp1^2}, amp1]       and               amplitude2[0] = 0.05; period2[0] = 4.010;     amplitude2[1] = 0.10; period2[1] = 4.015;     amplitude2[2] = 0.20; period2[2] = 4.02;     amplitude2[3] = 0.30; period2[3] = 4.035;     amplitude2[4] = 0.40; period2[4] = 4.05;     amplitude2[5] = 0.50; period2[5] = 4.075;     amplitude2[6] = 0.60; period2[6] = 4.105;     IntData = Table[{amplitude2[n], period2[n]}, {n, 0, 6}];     Fit[IntData, {1, amp2, amp2^2}, amp2]      and then tried this:               Plot[4.02 - 0.00698*amp2 + 0.271*amp2^2, {amp2, 0, 0.6},       AxesOrigin -> {0, 0},       2.84 - 0.00005*amp1 + 0.176*amp1^2, {amp1, 0, 0.6},       AxesLabel -> {"Amplitude", "Period"}, AxesOrigin -> {0, 0}]      thinking that would pu them both on the same graph, but I got an error message:               Plot::nonopt: Options expected (instead of {amp1,0,0.6}) beyond position 3 in Plot[4.02 -0.00698 amp2+0.271 amp2^2,{amp2,0,0.6},AxesOrigin->{0,0},2.84 -0.00005 amp1+0.176 amp1^2,{amp1,0,0.6},AxesLabel->{Amplitude,Period},AxesOrigin->{0,0}]. An option must be a rule or a list of rules. >>      Which I can sort-of-decipher but not really. So I am asking here what I did wrong. It's probably something obvious to many people here, but I am a rank beginner with Mathematica and am hoping that someone here can give me the kind of explanation they'd give to their grandmother. :-) No, I am not smarter than a 5th grader in this context. I hope this isn't too long a post. I'm just trying to get this done and offer my sincere thanks to those that answer. (And maybe this will help other newbies too, I can't be the only one)


title: What ways are there to convert an Integer to a list of digits?

text: How can I convert an `Integer` to list of numbers? For example:               781049 ==> {7, 8, 1, 0, 4, 9}      What are the possible ways of doing this? How well they compare performance- wise? Is a compiled version faster than `IntegerDigits`? How to convert huge amounts of large integers to digit form?


title: How to Open Kernel?

text: My MathKernel.exe resides in the following path: C:\Program Files\Wolfram Research\MathKernel.exe. I built a webMathematica application named "WebMathematica" and placed my entire application folder under the "Web-Apps" folder of Apache Geronimo. It's working fine. My Mathkernel was initialized using "mspConfiguration.xml" in the following way.               <KernelPool>       <KernelPoolName>General</KernelPoolName>         <KernelExecutable>C:\Program Files\Wolfram Research\MathKernel.exe</KernelExecutable>         <URLPattern>/*</URLPattern>     </KernelPool>      After this, I copied this entire folder and pasted it in I:\Projects. I opened the command prompt (cmd) and converted my application into a ".war" file. ![enter image description here](http://i.stack.imgur.com/iKYjj.jpg) Finally, I got a WebMathematica.war file. I then deployed this file in Apache Geronimo. I started my application server, and here I observed that Mathkernel wasn't initialized, so my application was not working as well. The problem was the kernel path, the one that I list in the "MSPCOnfiguration.xml" file. I believe that once a`.war` file has been created, we can't access outside paths. My question is how to fix this. I never use any `IDE`.


title: strange output from ToExpression

text: I want to set constraints to my parameters, e.g., "m", but Mma gives nonsensical output. How could I avoid this phenomenon by using another but "equivalent expression"? Or, is this just in my computer (i.e. a non- reproducible problem)? case 1. output is as expected:               ToExpression[      ToString[2.] <> ToString[">"] <> ToString["m"] <> ">=" <>        ToString[1.]]      > >     2. > m >= 1. >   case 2. output is unexpected:               ToExpression[      ToString[0.10465923419974825] <> ToString[">"] <> ToString["m"] <>        ">=" <> ToString[-6.0857416503171394*^-6]]      > >     -60.8574 >


title: Runge-Kutta-2 on System

text: After spending some time using the Mathematica documentation and this Mathematica.SE _answer_, I implemented the Runge-Kutta-2 routines. I am hoping someone can validate what I did and tell me that it is correct (especially the Butcher Tableau I used) and the step size $h = 0.1$.                 ClassicalRungeKuttaCoefficients[4, prec_] :=         With[{amat = {{1/2}}, bvec = {0, 1}, cvec = {1/2}},          N[{amat, bvec, cvec}, prec]]            {xf, yf} = {x, y} /.          First@NDSolve[{x'[t] == -y[t], y'[t] == x[t], x[0] == 1,           y[0] == 0}, {x, y}, {t, 0, 6},           Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4,          "Coefficients" -> ClassicalRungeKuttaCoefficients},           StartingStepSize -> 1/10];            xl = MapThread[Append, {xf["Grid"], xf["ValuesOnGrid"]}]            yl = MapThread[Append, {yf["Grid"], yf["ValuesOnGrid"]}]      We can find a closed form solution to this problem as:                 s = DSolve[{x'[t] == -y[t], y'[t] == x[t], x[0] == 1, y[0] == 0}, {x, y}, t]      Lastly, is there an automated way to update and step through each variant of RK-2, RK-3, RK-4, ... without having to manually enter the Butcher values (in other words, I want to step through each variant of RK-n and compare the errors (a table of that would be great))? There is a hint of this at _Wolfram's Reference page_. Note: I am embarrassed to say that I did not totally understand ClassicalRungeKuttaCoefficients (other than the coefficients).


title: Problem with NIntegrate when WorkingPrecision is specified

text: I am trying to evaluate this integral numerically: $$ \int_0^{\infty } m \exp (-m) J_1(m){}^2 \, dm $$ Everything is OK when only the integration method is specified:               NIntegrate[-m Exp[-m] BesselJ[1, m]^2, {m, 0, Infinity}, Method -> "ClenshawCurtisRule"]      but when I specify the `WorkingPrecision`, the integral remains unevaluated:               NIntegrate[-m Exp[-m] BesselJ[1, m]^2, {m, 0, Infinity}, Method -> "ClenshawCurtisRule",       WorkingPrecision -> 10]      What is wrong with this code? I am using _Mathematica_ v9.0.1 **UPDATE** This bug is still present in version 10.0.0.0.


title: How can I stop DiscreteMarkovProcess[] from relabeling vertices?

text: I'm attempting to calculate the mean first passage time between two vertices, $v_1$ and $v_2$, provided some undirected graph $G$. However, I noticed that running the following script:               DMPG = DiscreteMarkovProcess[1, G];     Print[Mean[FirstPassageTimeDistribution[DMPG,2]]      Returns a value for the mean first passage time between $v_1$ and some $v_i$ (instead of the desired $v_2$) due to a vertex name reassignment by DiscreteMarkovProcess[]. How can I stop this, or find the mapping between vertices in $G$ and nodes in the Markov chain generated by DiscreteMarkovProcess[]?


title: Alternative to using global variables in functions?

text: Often times I find myself using unevaluated variables (i.e. they show up blue) as dummy variables in whatever I'm returning such that I can evaluate them as needed later on. A simple example is as follows. Consider we have a function that returns a collection of monomials with 'x' as the independent variable.               getFns[n_] := Array[x^# &, n];      In this function, 'x' is simply a dummy (global) variable. Then suppose we have a function that returns a list of values corresponding to the list of monomials evaluated at specified coordinates:               useFns[fns_, coords_] := Array[fns[[#]] /. x -> coords &, Length[fns]];      Here is some sample input:               myFns = getFns[3]; (* {x, x^2, x^3} *)     myCoords = {1, 2, 3, 4};     useFns[myFns, myCoords]      which returns (as expected):               {{1, 2, 3, 4}, {1, 4, 9, 16}, {1, 8, 27, 64}}      The problem with this approach is that if the initial dummy variable is changed, it needs to be changed each time the list of functions is evaluated at specific points which can be tedious. Is there a better way to implement this?


title: Replacing a function in defined expressions

text: I'm having a little trouble with substituting functions into an expression. A minimal example is as follows:               diffy := D[U[z], z];     Unevaluated[diffy]/. U[z] :> Exp[-I*k*z]         Derivative[1][U][z]      Instead I want to do this:               Unevaluated[D[U[z], z]] /. U[z] :> Exp[-I*k*z]         -I E^(-I k z) k      but obviously I can't do the latter because my practical "diffy" expression is a result of many steps and consists of many lines of functions and derivatives of U[z]. Is there a way _Mathematica_ can take a function, substitute it into an expression such as "diffy" and evaluate it? I tried enclosing my expression above in "Evaluate[]" but that did not work. **Update:** The solution provided works for the above prototype, but when I add another level, it all breaks down. This is a minimum example:               U0 := U[z];     diffy1 := Hold[D[U0, z]];     diffy1 /. U[z] :> Exp[-I*k*z] // ReleaseHold         Derivative[1][U][z]      I have tried holding the first expression, releasing the hold in the definition of the second one and then holding again - and many other combinations of holding/releasing but nothing has worked so far. Is there a fundamental reason why this is not working? Maybe I'm not thinking of Mathematica correctly but I don't understand why this does not work immediately.


title: Select elements in a list when the firsts two items are the same (no matter order) and its third element is lower

text: I'm trying to solve a heuristic search in Mathematica 9. I get a list of succesor nodes called open with two nodes and a heuristic value, e.g.,               open={{node1,node2,H1},...,{node2,node1,Hi},...,{noden,noden,Hn}}      I need to select only when the two firts elements are the same (no matter order) but where the heuristic is lower. Following the code that I used for sort the elements by the first item               Sort[open, #1[[1]] < #2[[1]] &]      I try to select cases when the same first and second numbers and the heuristic value where lower:               Select[open, #1[[1]]==#2[[1]] &&  #1[[2]] == #2[[2]] && #1[[3]]< #2[[3]] &]      But, returns `{}`... My list               open = {{1, 2, 3}, {2, 1, 5}, {1, 2, 2}, {3, 4, 5}, {4, 3, 2}};      Between `{1, 2, 3}`, `{1, 2, 2}`, and `{2, 1, 5}` then choose `{1, 2, 2}` that has the lower heuristic. And, between `{3, 4, 5}` and `{4, 3, 2}` then `{4, 3, 2}`...               {{1, 2, 2},{4, 3, 2}}      Thanks,


title: Accessing Draw Tools objects with notebook commands

text: In a notebook I would like to be able to access and modify the properties of draw tools objects after I've drawn them on a graph by hand. For example, if make a plot and, with draw tools, make a circle, is it possible now to access the circle Graphics object and its properties with commands in the notebook?


title: How to create a 'shortest-path link matrix'?

text: I am currently struggling with a problem in _Mathematica_. I have a randomly generated weighted Barabasi-Albert graph and a corresponding demand matrix of the same size. I'm generating the graph with the following code:               generateBarabasiAlbertGraph[n_, l_] /; n >= 3 :=       Module[{g = CycleGraph[3], vc, vl, el},      Do[      vc = VertexCount[g];      vl = VertexList[g];      el = EdgeList[g];      ec = EdgeCount[g];      wl = RandomInteger[{5, 20}, {ec}];      graph = Graph[el, EdgeWeight -> wl, EdgeWeight -> "EdgeWeight"];      el = Map[UndirectedEdge[vc + 1, #] &,         RandomSample[VertexDegree[g] -> vl, Min[l, vc]]        ];      g = Graph[Join[EdgeList[g], el],         VertexLabels -> Placed["Name", Before]], {n - 3}];        g;      ];     v = RandomInteger[{5, 10}];     l = RandomInteger[{1, 2}];     generateBarabasiAlbertGraph[v, l];      I then calculate the shortest path from each vertex `a` to vertex `b` if vertex `a` has a demand from vertex `b`. Now I would like to create a _shortest-path link matrix_ for all shortest paths in the graph. (A similar graph appears in this paper by Kelly). This is what I already figured out:               TableForm[Array[0 &, {EdgeCount[graph], 1}], TableHeadings -> {EdgeList[graph], {sp}}]      ![enter image description here](http://i.stack.imgur.com/Eh87P.png) As you can see, I have the shortest-path 3-1-4 and need to replace the 0 at the corresponding row if this edge is used in the shortest path. So in this case I need a `1` in row 3 and 4 Any ideas? Best, Julian


title: 3D Plot from 3 Polynomial Equations

text: Here is my code:               f0[\[Lambda]0_, y0_, y1_, x_, k_] :=-(-y0 + k*x*y1 + (1 - k*x^2)*\[Lambda]0^2)^2 + 4*x^2*(1 - k^2*x)*\[Lambda]0^2*Abs[-(k*y0) + y1 + k*\[Lambda]0^2 - x*\[Lambda]0^2]          f1[\[Lambda]1_, y0_, y1_, x_, k_] := -(x*y0 - y1 + (1 - k*x^2)*\[Lambda]1^2)^2 + 4*x^2*(-k + x)*\[Lambda]1^2*Abs[-(k*y0) + y1 - \[Lambda]1^2 + k^2*x*\[Lambda]1^2]          f2[\[Lambda]0_, \[Lambda]1_, x_, k_] :=-1/27 + \[Lambda]0^3 + 3*x*\[Lambda]0^2*\[Lambda]1 + 3*k*x*\[Lambda]0*\[Lambda]1^2 + \[Lambda]1^3      $k$ is some positive number which can be chosen freely for example $k=1$ or $k=3/2$ or $k=2/3$. $f_0$ is only a function of $\lambda_0$, likewise $f_1$ is only a function of $\lambda_1$. If I would insert $\lambda_0$ from $f_0$ into $f_2$ and insert $\lambda_1$ from $f_1$ into $f_2$ then $f_2$ will be a function of only $x,k,y_0,y_1$. In this case I would like to have a 3D plot with x axis =$y_0$, y-axis=$y_1$, and z-axis =$x$ for $7/81<y_0<1/9$ and $7/81<y_1<1/9$, for some known $k$ as mentioned above and $x$ must be eventually in $0<x<1$, by positivity condition. Is it doable? how can I do that? **EDIT:** I obtain the polynomials via equating them to $0$. Namely, $f_0=0$, $f_1=0$ and $f_2=0$ are all given and known. Thank you very much.


title: Scripting a visual comparison between the derivative and some finite scheme

text: I want to make an application where I can specify a function $f:\mathbb R\to \mathbb R$ and a dynamic value $n$ and get the function $\frac{f(x+\frac{1}{n})-f(x)}{\frac{1}{n}}$. I moreover want to compare it to $f'(x)$, i.e. plot it in the same graph. I first tried it with $f(x):=x^m, m=2$ and used Listplot, but even here I got into some obsticles. Here is my attempt:               Table[((x+1/n)^m)-x^m)/(1/n),{m,1000,100}]          Function[x,#][Rangle[10]]&/@%          %/.{n->3}//N          Show[     Plot[3 x^2,{x,0,10}]     ListPlot[#]&/@%     ]      If I do just _ListPlot[#] &/@%_ instead of the last line, the I get all the plots, however it doesn't compare. I'd also like to get away from $x^m$ and make an abstraction w.r.t. $f$, but my tries with Function fails. Is there a difference between _Function[$x$,term]_ and _$\lambda x.$term_ ?.


title: making a function of fit results

text: I have 13*25 datasets consisting each of around 40 points. They are all linear I can easily fit them using fit function but my problem is when I want to call them I want them to be functions. I want to define a function using fit results something like               Radius[s1,s2,o,i]=Fit[data[s1,s2,o],{1,x},x]      Do you have any idea how can I make this function? any ideas? Edit: My data is a set data of radius, {i,r(i)} what are s1,s2 and o, Here it come in details if you are interested ![Molecules](http://i.stack.imgur.com/GEjtR.jpg) s1 is site one on molecule one and s2 is site2 on molecule 2 an o is orientation meaning the way molecules can approach each other anyways in each case we will have different radiuses between these molecules and data[s1,s2,o] is a set of 40 points for example the picture is orientation one, data[1,7,1] will be 40 points of distance between atom 1 in molecule 1 and atom 7 in molecule 2.               {{1, 6.44479}, {2, 6.34609}, {3, 6.24742}, {4, 6.1488}, {5,     6.05023}, {6, 5.9517}, {7, 5.85322}, {8, 5.7548}, {9, 5.65642}, {10,     5.55811}, {11, 5.45986}, {12, 5.36167}, {13, 5.26356}, {14,     5.16551}, {15, 5.06754}, {16, 4.96965}, {17, 4.87185}, {18,     4.77414}, {19, 4.67652}, {20, 4.57901}, {21, 4.48161}, {22,     4.38432}, {23, 4.28716}, {24, 4.19013}, {25, 4.09325}, {26,     3.99652}, {27, 3.89995}, {28, 3.80356}, {29, 3.70736}, {30,     3.61137}}      the first iterator is the step number, meaning by each step molecules get closer to each other. Thanks for your attention


title: How to prevent Mathematica combining similar terms?

text: For example, given: `x + 2 x y - x y + y` I need to put it in this form: `x (1 + 2 y) - y (x - 1)`. But _Mathematica_ will simplify (automatically) `x + 2 x y - x y + y -> x + x y + y` and I will not be able to group this expression as I need (I need to factor `(1 + 2 y)`). The main question is how to prevent _Mathematica_ simplifying, _or_ how to factor out `x (1 + 2 y)` from this expression: `x + x y + y`? Alternatively, how to put `x + x y + y` into `x (1 + 2 y) - y (x - 1)`?


title: Computing the average over the values of the boxes surrounding a box in a grid

text: Recenty,I encounter a problem described as below: > I have a table that owns many boxes,and each box has a sequence like > 1,2,3,4,etc.In addition,every box has a value. > > Now I want to make the value of a box become the average of the values of > boxes that arounding it. > > For example, > > `the value of box 1` become the average of`the values of box 2,11,12`, > > `the value of box 11` become the average of`the values of box 1,2,12,22,21` > > `the value of box 14` become the average of`the values of box > 3,4,5,13,15,23,24,25`               Grid[Reverse@Partition[Range[80], 10], Frame -> All]      ![enter image description here](http://i.stack.imgur.com/ufgS6.png)               sequenceValue = Thread@List[Range[80], RandomInteger[{0, 1}, 80]];     Reverse@Partition[sequenceValue, 10] // Grid[#, Frame -> All] &      ![enter image description here](http://i.stack.imgur.com/Mizh4.png) ### My trial **Extract the sequence of the boundary of table**               ExtractBoundarySequence[length_, width_] := Block[       {bottomData, topData, leftData, rightData, cornerData},      bottomData = Range[2, length - 1];      topData = Range[2 + (width - 1) length, width*length - 1];      leftData = Range[1 + length, 1 + (width - 2) length, length];      rightData = Range[2 length, (width - 1) length, length];       cornerData = {1, length, 1 + (width - 1) length, width*length};      List[cornerData, bottomData, topData, leftData, rightData]     ]      **Extract the sequence of the innerior of table**               ExtractInteriorSequence[length_, width_] := Block[     {},     DeleteCases[      Flatten@Table[{# - 1, #, # + 1} + i *length, {i, -1, 1}], #] & /@      Flatten@Table[Range[2 + length, 2 length - 1]+(i - 1)length , {i, 1, width - 2}]     ]      **Caculate the value of the boundary of table**                CalculateBoundaryValue[SequenceValue_, length_, width_] := Block[      {corner, bottom, top, left, right,       cornerSequence, bottomSequence, topSequence, leftSequence, rightSequence,      value1, value2},      {corner, bottom, top, left, right} =ExtractBoundarySequence[length, width];      bottomSequence = DeleteCases[       Flatten@Table[{# - 1, #, # + 1} + i *length, {i, 0, 1}], #] & /@bottom;      topSequence = DeleteCases[       Flatten@Table[{# - 1, #, # + 1} + i *length, {i, -1, 0}], #] & /@top;      leftSequence = DeleteCases[       Flatten@Table[{#, # + 1} + i *length, {i, -1, 1}], #] & /@ left;      rightSequence = DeleteCases[       Flatten@Table[{# - 1, #} + i *length, {i, -1, 1}], #] & /@ right;          value1 =      Mean /@ (Flatten[       Map[Part[SequenceValue[[All, 2]], #] &,        {bottomSequence, topSequence, leftSequence,          rightSequence}, {2}], 1]);          cornerSequence =     {      DeleteCases[        Flatten@Table[{#, # + 1} + i*length, {i, 0, 1}], #] &[corner[[1]]],      DeleteCases[        Flatten@Table[{# - 1, #} + i *length, {i, 0, 1}], #] &[corner[[2]]],      DeleteCases[        Flatten@Table[{#, # + 1} + i *length, {i, -1, 0}], #] &[corner[[3]]],      DeleteCases[        Flatten@Table[{# - 1, #} + i *length, {i, -1, 0}], #] &[corner[[4]]]     };          value2 =       Mean /@ (Part[SequenceValue[[All, 2]], #] & /@ cornerSequence);           Thread@List[Flatten@{corner, bottom, top, left, right}, Join[value2, value1]]     ]      **Caculate the value of the innerior of table**                CalculateInteriorValue[SequenceValue_, length_, width_] := Block[      {interiorSequence, interiorNumber, value},      interiorSequence =      ExtractInteriorSequence[length, width];      interiorNumber =       Flatten@Table[Range[2 + length, 2 length - 1]+(i - 1)length, {i, 1, width - 2}];      value =          Mean /@ ((Part[SequenceValue[[All, 2]], #] &) /@ interiorSequence);       Thread@List[interiorNumber, value]      ]      **Final sequence-v-alue**               NewSequenceValue[SequenceValue_, length_, width_] := Block[      {NewBoundaryData, NewInteriorData},      NewBoundaryData =      CalculateBoundaryValue[SequenceValue, length, width];      NewInteriorData =      CalculateInteriorValue[SequenceValue, length, width];      SortBy[Join[NewBoundaryData, NewInteriorData], First]     ]      **Using my function:**               Reverse@Partition[NewSequenceValue[sequenceValue, 10, 8], 10] //Grid[#, Frame -> All] &      ![enter image description here](http://i.stack.imgur.com/GyJOm.png) and it can achieve the result. However, I think my method is a little tedious, so my question: is there a better algorithm (method) to solve my problem?


title: Add noise (random values) to data

text: I have a data set(x,y,z) corresponding to a hemisphere, to which I want to add some controlled noise (random values) to make the x,y,z values diverge. > data = Flatten[ Table[{i, j, i^2 + j^2}, {i, -10, 10, 1}, {j, -10, 10, 1}], > 1] The above is a 3 coordinate set (x,y,z) values for a hemisphere. How do I add random noise to it? Thank you


title: Classify command in Mathematica 10

text: Where can I find documentation about the methods used in the `Classify` command in _Mathematica 10_? Specifically, how does the method `Markov` work? Any pointers to the documentation for the methods available in `Classify`, or an explanation of how the method `Markov` works, will be appreciated.


title: Simulating a bounded continuous random walk in $n$-dimensions

text: You're given three $n$-dimensional real vectors, $\mathbf{x_0}$, $\mathbf{a}$ and $\mathbf{b}$, with $\mathbf{a} \le \mathbf{x_0} \le \mathbf{b}$ (vector inequality means component-wise inequality). Starting at $\mathbf{x_0}$, simulate a bounded continuous random walk, such that the end point $\mathbf{x}$ of the walk is always inside the bounds $\mathbf{a} \le \mathbf{x} \le \mathbf{b}$. I want to program this with Mathematica. The parameters of the simulation are the bounds $\mathbf{a}, \mathbf{b}$, the starting point $\mathbf{x_0}$, and a parameter $\sigma$ that controls the standard deviation of the step size. **What's a _simple_ and _efficient_ way to do this?** **Here's what I've tried:** The following function `Walk` takes a vector $\mathbf{x}$ and returns the next point on the random walk.               Walk[x_List, a_List, b_List, sigma_] :=       x + RandomVariate[        TruncatedDistribution[Transpose[{a - x, b - x}],         MultinormalDistribution[ConstantArray[0, Length[x]],           sigma*IdentityMatrix[Length[x]]]]]      For some reason this is exceedingly slow. You can try it yourself. After executing:               a = ConstantArray[0, 10];     b = ConstantArray[1, 10];     x = RandomReal[{0, 1}, 10];      The following line of code takes forever on my PC:               Walk[x, a, b, 4]      **Why?**


title: Stacked bars in a barchart with individual and overall labels

text: There are a lot of questions here about formatting of barcharts, but I couldn't find one matching mine. I seek to display data in stacked format, labeling each component of the stack and the total for each bar. I can't make it work for the life of me. As a starting place, consider               BarChart[{{1, 2}, {2, 3}, {2, 4}}, ChartLayout -> "Stacked",       ChartLabels -> {Range[2004, 2006], None}]      I want each component labeled as below -- ![goal](http://i.stack.imgur.com/Wr6C9.png) My efforts to date result in some pieces appearing but never all, and often in the wrong places.


title: Use output of NDSolve as an integration limit of NIntegrate

text: I've been working on a ODE, for a function `y[t]`, that, due to it's complexity, has to be solved numerically. I then need to use the results of NDSolve, which I store in the variable $s$, in the following manner:               Plot[          NIntegrate[                     f[r],{r,b,Evaluate[Sqrt[{(y[t])^2+b^2}/.s]]}                    ],          {t,0,5}         ]      which returns the following errors:               NIntegrate::nlim: r = {20.025} is not a valid limit of integration. >>     NIntegrate::nlim: r = {20.025} is not a valid limit of integration. >>     NIntegrate::nlim: r = {19.6024} is not a valid limit of integration. >>     General::stop: Further output of NIntegrate::nlim will be suppressed during this calculation. >>      It seems that the problem is due to $\text{Evaluate}$ returning the number in parentheses...I've tried the function $\text{N}$ but it has the same problem. Maybe I'm doing something else wrong. Can anyone help me? **Edit:** I shall post the complete code as suggested in the comments:               Module[{b = 1., k = 1., Vt0 = 2., L0 = 20.},      s = NDSolve[        {y''[t] - (y'[t])^2 (b^2/((b^2 + (y[t])^2) y[t]) - y[t]/(              b^2 (Exp[(y[t])^2/b^2] - 1))) - (z'[t])^2 (-1/(              y[t] Exp[(y[t])^2/b^2]) (Exp[(y[t])^2/b^2] -                 1) (b^2 + (y[t])^2)) + k*(b^2 y[t])/(b^2 + (y[t])^2) == 0,          y[0] == L0, y'[0] == 0,         z''[t] + 2*y'[t]*z'[t]*(y[t]/(b^2 + (y[t])^2)) == 0, z'[0] == Vt0,          z[0] == 0},        {z, y}, {t, 0, 5}        ];      f[r_] := Sqrt[E^(r^2/b^2)/(-E + E^(r^2/b^2)) - 1];           Plot[       NIntegrate[f[r], {r, b, Evaluate[Sqrt[(y[t])^2 + b^2] /. s]}],       {t, 0, 5}       ]      ]


title: How to use Button as a control in Manipulate

text: I am trying to use `Button` as a controller in `Manipulate`. I have done some thing like this:               Manipulate[       p = Column[{Plot[Sin[x - n], {x, 0, 12.4}, ImageSize -> 500]}],        {n, 0, 1},        Button["Export", Export["test.jpg", p]]]      Is there any better way to do this? I am thinking of making something like this:               Manipulate[       Column[{p = Plot[Sin[x - n], {x, 0, 12.4}, ImageSize -> 500],        If[x, Export["test.jpg", p]]}],        {n, 0, 1},        {x, {False, True}}]      in which the controller `x` is represented by `Button` instead of a `Checkbox`. How can I do that?


title: WSTP Callbacks into CDF

text: Essentially, what I am trying to find out is the limitation of the WSTP/LINK protocols in Enterprise-licensed CDF. Specifically, are any callbacks to the kernel from those protocols allowed in CDF , Enterprise CDF, or Player Pro? Are any methods disabled? This problem is difficult since very few people discuss the issue, and Wolfram documentation makes little mention of disabling specific functionality for CDF/Player Pro (e.g. http://www.wolfram.com/player-pro/how-player-pro- compares.html is not very specific). More importantly, if callbacks were allowed without any restrictions on what `Evaluate` could do, what would stop a person from evaluating arbitrary expressions using a combination of a LINK and IPC (i.e. interprocess communication)? I have been told by Wolfram representatives that version 10 Enterprise CDF should have less LINK restrictions that the previous version (e.g. .NET/LINK is allowed now), but that still doesn't help me understand where Wolfram draws the line. I appreciate any answers or suggestions. If I hear back from Wolfram on the issue, I will update this post.


title: Design considerations behind `O` (a.k.a. BigOh, a.k.a. Landau Order)

text: This works without any warnings: `O[Log[x]]`. This raises a warning: `O[x^2]`. I have a few questions around this:   1. Why is it a good idea for `O[x^2]` to raise a warning when `O[Log[x]]` does not?   2. Is `O[Log[x]]` meaningfully useful in _Mathematica_ given the semantics of its `O[]`?   3. How does one represent $\mathcal{O}(\log(x-x_o))$ naturally in _Mathematica_? (Probably there is no way because it is useless...)   4. Why is `O[x, x0]^2` the representation chosen by _Mathematica_ authors instead of something like `O[(x-x0)^2]`?


title: How to neglect higher power terms in a polynomial expression

text: I have a polynomial expression of order n (say n=20). `F(x)=1+x+x^2+x^3++...x^20`. I want to approximate the polynomial for order 3 only. So I need to make the coefficient of power of x greater than 3 to zero. How can I do that in _Mathematica_?


title: Computing the fundamental matrix and its monodromy matrix

text: Assume $\dot x(t)=v(x,t)$ is a $T$-periodic, with respect to $t$ dynamical system. That is $x\in\mathbb R^n,\quad t\in\mathbb R,\qquad v(x,t+T)=v(x,t)$. Let $x_0$ be a periodic solution, and in elementary functions. Assume if you must that $x_0$ is as smooth as you need. How can one use Mathematica to compute symbolically a fundamental matrix for the linearised around $x_0$ system and its monodromy matrix? Any help would be appreciated.


title: How can I monitor and detect when a website has changed?

text: I'd like to detect when a website has changed, in particular when new products are announced on a site that doesn't change to often. Here's the code I wrote but I don't think it is perfect. Also the site might block me if I ping it too many times...               content = Hash@URLFetch["www.wolfram.com"];     $c = 0; Dynamic[$c]          RunScheduledTask[$c++;       If[content =!= Hash@URLFetch["www.wolfram.com"],        SendMail["To" -> "me@me.com",         "Subject" -> "The wolfram site has been updated!",         "Body" -> {"Check wolfram now"},         "From" -> "email", "Server" -> "smtp.gmail.com",         "UserName" -> "email", "Password" -> "pass",         "PortNumber" -> 587, "EncryptionProtocol" -> "StartTLS"];          content = Hash@URLFetch["www.wolfram.com"]], 100]      I'd also like to create a standalone executable for running this process in the background without mathematica being open.


title: Using MathKernel without echoing a window

text: I use M9 under windows 7 on a command prompt as: MathKernel -noprompt -script file.m On execution, there is a like a kernel window opened and then closed immediately. How can I avoid this window and run my script in silence? thanks


title: Creating a simple function to compute the average of the difference between pairs of elements in an array

text: I have a list of the form:               {{{1, 2}, {3, 4}}, {{5, 6}, {7, 8}}, {{9, 10}, {11, 12}}}      How can I create a function which will take each set at the highest level of the list:               {{1, 2}, {3, 4}}          {{5, 6}, {7, 8}}          {{9, 10}, {11, 12}}      And subtract the first element of the set from the second element of the set, here:               {3, 4} - {1, 2} = {2, 2}          {7, 8} - {5, 6} = {2, 2}          {11, 12} - {9, 10} = {2, 2}      Then compute the average or median of these values (`{2, 2}` here in either case)? It's clear how to do this using loops or tables, but is there a simpler way to proceed that makes clever use of _Mathematica_? For a table based implementation:               PairList = {{{1, 2}, {3, 4}}, {{5, 6}, {7, 8}}, {{9, 10}, {11, 12}}};     PLMean = Mean[Table[PairList[[u, 2]] - PairList[[u, 1]], {u, 1, Length[PairList]}]]


title: Rotate lines around their midpoint

text: I am trying to rotate a series of lines by the same angle around their centre point. Say each line is defined by a co-ordinate at each end, simplified example below of a typical dataset.               foo = {{{100, 25}, {150, 45}}, {{200, 45}, {240, 85}}};      Then the midpoint can be found by:-               boo = (#[[2]] - #[[1]])/2 + #[[1]] & /@ foo;      This is visualised using the code below:               Show[     Graphics[{Thick, Dashed, Black, Line /@ foo}],     Graphics[{Black, PointSize[Large],     Point[boo],     Table[{Black, Text[ToString[boo[[p]]], boo[[p]] + 2]}, {p, Length[boo]}]}]     ]      ![enter image description here](http://i.stack.imgur.com/zmjPE.jpg) I was hoping to create a series of transformations using RotationTransform such as below:               rot = RotationTransform[Pi/2, #] & /@ boo      Then mapping this across foo using MapThread, but I'm not sure which function I should be mapping across to make each line rotate by 90 degrees around its own centre point.               MapThread[?, {woo, rot}]      I think I'm close to what I want to achieve (or maybe not), but I can't quite seem to close it out. All suggestions welcome.


title: Once more on object orientation in Mathematica: does it have to be so hard?

text: Recently I came across a set of problems which would be solved most easily within an object-oriented approach. I first attempted to solve them by other means, but found the complexity of the code growing too fast and, while individual functions used a mixture of styles, the overall feel and organization of the code was de-facto procedural. Having looked at the final code, I saw that I would have written it much faster and cleaner in Java, which is when I realized that OO was needed. Basically, I needed it to reuse some behavior and decrease coupling between components. We all know that there were a number of attempts on implementing OO extensions for Mathematica. However, none of them IMO were very simple, and none of those I looked at had a natural "look and feel" for Mathematica (I could have missed some). OTOH, we also know that Mathematica has powerful metaprogramming capabilities, which should make it possible to implement the core of an OO extension rather easily. Since OO is a broad paradigm, here are some requirements to restrict the problem. The desired OO extension should   * _Be idiomatic_. Programming in it should feel natural for experienced Mathematica users. In particular, it should go well with immutable Mathematica expressions, and have the least possible problems concerning interoperability with Mathematica (garbage collection - related problems etc).   * _Support instantiation, inheritance and polymorphism_. By inheritance I mean reuse of behavior (methods) rather than state (fields).   * _Reuse_ as much of the Mathematica's core constructs as possible   * _Limit the things that can be done with it as little as possible_ , compared to the programming techniques and tricks we are used to in Mathematica   * _If possible, have at least minimally convenient syntax_ for things like method calls.   * _Implementation should be simple_ in the sense that it should not introduce too many exterior constructs, such as `ToString` \- `ToExpression` cycles, needless manipulations with contexts, creation of new symbols which are not strictly necessary, etc.  I am not so much worried about efficiency here, since the intended purpose for such an extension is more to help with larger-scale code organization for certain classes of tasks, than to be used with millions of light-weight objects. In many previous questions on similar topics, the two rather separate topics were often mixed together. One concerns mutable struct-like data structures and their possible implementations in Mathematica, and the other was about OO proper. So, to clarify a bit more, I am not so much interested in the former topic here. I am rather more interested in the dynamic aspects (behavior reuse), and ways to conveniently package code and organize larger projects, but not so much in efficient mutable data structures. In other words, the subject matter is what can the object orientation bring to the table for larger Mathematica projects, and how can we make the most out of it in Mathematica. **So, the question is** : can we have a core of such an OO extension in under 100 lines of Mathematica code?


title: How do I apply a rule to every element in a matrix?

text: I need to apply a set of rules (I used `Which`) to every element in a matrix, but the rules are dependent on the position `[[i, j]]` of the matrix.


title: Improved interpolation of mostly-structured 3d data

text: This question arose in response to a comment by Leonid to my answer for this question. He noted that for unstructured grids, `Interpolation` can only use `InterpolationOrder->1`. For example:               data = Flatten[Table[{x, y, x^2 + y^2}, {x, -10, 10}, {y, -10, 10}], 1];     dataDelete = Delete[data, RandomInteger[{1, Length[data]}]]     intD = Interpolation[dataDelete]     (* Interpolation::udeg: Interpolation on unstructured grids is currently only supported for InterpolationOrder->1 or InterpolationOrder->All. Order will be reduced to 1. *)      which gives a result much worse than `All`-order `Interpolation`. So, here's my question: For data that is _largely_ structured, but is missing one (or a few?) grid points, is there any way to figure out where the grid is missing points? This would allow you to first use linear interpolation to find a good guess for the value of the function at the missing grid points, then add that grid point to the dataset, only to interpolate the whole thing after. Then the error in the interpolation will be localized to the reconstructed region. Something like               dataReconstructed=Append[dataDelete,Sequence @@ Flatten@{#, intD @@ #} & /@missingcoords]     intReconstructed = Interpolation[dataReconstructed]      A manual example:               todelete = RandomInteger[{1, Length[data]}];     dataDelete = Delete[data, todelete];     intD = Interpolation[dataDelete]     missingcoords = {data[[todelete, {1, 2}]]}          dataReconstructed = Append[dataDelete,Sequence @@ Flatten@{#, intD @@ #} & /@missingcoords]     intReconstructed = Interpolation[dataReconstructed]      Comparing these two methods using `exact[x_, y_] := x^2 + y^2`:               Plot3D[intD[x, y] - exact[x, y], {x, -10, 10}, {y, -10, 10}]      ![Mathematica graphics](http://i.stack.imgur.com/WJAgI.png) has differences all over the interpolated region, worse at the missing point. But:               Plot3D[intReconstructed[x, y] - exact[x, y], {x, -10, 10}, {y, -10, 10}]      ![Mathematica graphics](http://i.stack.imgur.com/n4pSD.png) is much better everywhere except for the missing points. In order to do this for a set where I don't know which grid points are missing, **is there a way to figure out where the missing grid points are in a mostly structured grid?**


title: Changing color of an object in an image

text: Consider the following image: ![Red ball in a non-red environment ](http://i.stack.imgur.com/Qr7Tx.jpg) How can I change all red colors in this image into (for example) blue.


title: Small Issue with Chebyshev Derivative Appoximation

text: I am trying to get approximate the derivative of a function from its Chebyshev expansion. I start out with the following random function               In[4]:= f[x_] = -4.1261 + .99213*x - 0.562172*x^2 - 3.7120 x^3 + 3.06508*x^4      Which looks like ![graph of random function](http://i.stack.imgur.com/v4aNt.jpg) Now I try the Chebyshev approximations as follows:               In[105]:= cgs = 5          In[107]:= Clear[A, a]          In[108]:= A = Array[a, cgs + 1, {0}]          Out[108]= {a[0], a[1], a[2], a[3], a[4], a[5]}          In[109]:=   a[0] = (1/Pi)* Integrate[(f[x]* ChebyshevT[0, x]/Sqrt[1-        x^2]), {x, -1, 1}]          Out[109]= -3.25778          In[110]:= For[i = 1, i <= cgs, i++,  a[i] = (2/Pi)*Integrate[(f[x]*      ChebyshevT[i, x]/Sqrt[1 - x^2]), {x, -1, 1}]]          In[111]:= A          Out[111]= {-3.25778, -1.79187, 1.25145, -0.928, 0.383135, -1.97901*10^-15}          In[113]:= cf = Table[a[i]*ChebyshevT[i, x], {i, 0, cgs}]          Out[113]= {-3.25778, -1.79187 x, 1.25145 (-1 + 2 x^2), -0.928 (-3 x + 4          x^3), 0.383135 (1 - 8 x^2 + 8 x^4), -1.97901*10^-15 (5 x - 20 x^3 + 16 x^5)}          In[115]:= CF[x_] = Sum[cf[[i]], {i, 1, cgs + 1}]          Out[115]= -3.25778 - 1.79187 x + 1.25145 (-1 + 2 x^2) - 0.928 (-3 x + 4 x^3)          + 0.383135 (1 - 8 x^2 + 8 x^4) -  1.97901*10^-15 (5 x - 20 x^3 + 16 x^5)          Plot[{CF[x], f[x]}, {x, -1, 1}]      gives ![graph of function and Chebyshev expansion](http://i.stack.imgur.com/p3m98.jpg) So I think I have done everything correctly up until now, because I cannot see a difference between the two graphs. As for the derivative approximation, I use the discussion in Numerical Recipes section 5.9. First I take the derivative of `f[x]`               In[200]:= df[x_] = D[f[x], x]          Out[200]= 0.99213 - 1.12434 x - 11.136 x^2 + 12.2603 x^3          In[201]:= A          Out[201]= {-3.25778, -1.79187, 1.25145, -0.928, 0.383135, -1.97901*10^-15}          In[236]:= Clear[dA, da]          In[237]:= dA = Array[da, cgs + 1, {0}]          Out[237]= {da[0], da[1], da[2], da[3], da[4], da[5]}      Now I use equation (5.9.2)               In[273]:= da[cgs] = da[cgs - 1] = 0;          In[275]:= For[i = 4, i >= 1, i--, da[i - 1] = 2*(i)*a[i] + da[i + 1]]      to get               In[276]:= dA          Out[276]= {-9.15174, 8.0709, -5.568, 3.06508, 0, 0}          In[277]:= dcf = Table[da[i]*ChebyshevT[i, x], {i, 0, cgs}]          Out[277]= {-9.15174, 8.0709 x, -5.568 (-1 + 2 x^2), 3.06508 (-3 x + 4 x^3),            0, 0}          In[278]:= dCF[x_] = Sum[dcf[[i]], {i, 1, cgs + 1}]          Out[278]= -9.15174 + 8.0709 x - 5.568 (-1 + 2 x^2) + 3.06508 (-3 x + 4 x^3)          Plot[{dCF[x], df[x]}, {x, -1, 1}]      which gives ![plot of wrong derivative](http://i.stack.imgur.com/vmRU4.jpg) which is clearly not the same graph. I was able to figure out that the problem is in the constant term. Namely, it is twice the value it should be. To test this idea I created the `test` function               In[279]:=  test[x_] = -9.151740000000007`/2 + 8.070896000000015` x -      5.568000000000007` (-1 + 2 x^2) + 3.065080000000011` (-3 x + 4 x^3);      and plotted `Plot[{test[x], df[x]}, {x, -1, 1}]` to get ![plot of derivative](http://i.stack.imgur.com/64XHi.jpg) which shows that `test` is the same as `f[x]`. If anyone could help spot my mistake, that would be greatly appreciated.


title: Return a Part by Reference? (For accessing hierarchical / structured data)

text: Is there a way to return a part of a list by reference, for reading and writing? I want to do this so I can easily manipulate "rule-styled structured", i.e., hierarchical data that is loaded from a MAT file with `struct`'s. Example:               obj = {"Name" -> "Something", "Age" -> 10};      I think this is a relatively simple way to get / set fields, which deals with nested structures by limiting `Position` to `levelspec` 1:               RuleQ[x_] := False;     RuleQ[rule_Rule] := True;     RuleQ[{}] := True; (* An empty list can be a Rule set... *)     RuleQ[rules_List] := SameQ[Head[First[rules]], Rule];          StructQ[x_] := False;     StructQ[rules_List?RuleQ] := True;          FindFirst[list_List, item_] := Module[{res},        res = Flatten[Position[list, item, 1]];        Return[If[Length[res] > 0, First[res], {}]];        ];          FieldRule[field_] := field -> _;          GetFieldIndex[obj_?StructQ, field_] := FindFirst[obj, FieldRule[field]];          SetField[obj_?StructQ, field_, value_] := Module[{index},        index = GetFieldIndex[obj, field];        If[NumberQ[index],         obj[[index, 2]] = value;         ,         AppendTo[obj, field -> value];         ];        ];     SetAttributes[SetField, HoldFirst];          GetField[obj_?StructQ, field_] := Module[{index},        index = GetFieldIndex[obj, field];        (* How to return a reference? *)        Return[obj[[index, 2]]];        ];     SetAttributes[GetField, HoldFirst];      So you can do something like               Print[GetField[obj, "Name"]];     SetField[obj, "Name", "REPLACED"];     Print[GetField[obj, "Name"]];      All that being said, I would like to do something like this:               GetField[obj, "Name"] = "REPLACED";      An example of my hope for such a feature:               GetFieldByReference[obj_?StructQ, field_] :=         Return[Part[obj, {GetFieldIndex[obj, field], 2}];     SetAttributes[GetFieldByReference, HoldReturn];     obj = {"Name" -> "Something", "Age" -> 10};     GetFieldByReference[obj, "Name"] = "REPLACED";      The contents of `obj` would then be:               {"Name" -> "REPLACED", "Age" -> 10}


title: Numerically obtaining the inverse Laplace transform of data

text: I have been using several Mathematica packages to do numerical inverse Laplace transforms on known (expressible in closed form) expressions, $\tilde{f}(s)$. I am now being confronted with the more difficult problem where I can only obtain $\tilde{f}(s)$ as a table like               {{s_1,f_1},{s_2, f_2}, {s_3, f_3 ... etc.}      I have not been able to find any packages for doing the inverse Laplace transform on tabulated data. I saw one ancient discussion in a forum that says the inverse Fourier transform can be used but I don't see how that's possible. I would greatly appreciate if anyone had done something like this before, or knew of an appropriate package, or had any suggestions for a workaround. EDIT: this is how I'm getting tabulated results for the transformed function: I have $\tilde{f}(s) =\frac{1}{s} \frac{1-\tilde{g}(s)}{1-a \tilde{g}(s)}$ where $\tilde{g}(s) = \int_{-\infty}^0 h(s,x) dx$. I know what $h(s,x)$ is, but the integral over $x$ cannot be done analytically so I am doing it numerically which is why I end up with tabular data for $\tilde{f}(s)$.


title: interpolation of 3D data

text: I have a 3D surface given in data-points of the form ${x,y,z}$. What is the easiest way to get the interpolated value $z=f(X,Y)$ for given coordinates ${X,Y}$ (which are of course not in the data list)?


title: Calculating mean of multiple stochastic processes and setting a lower threshold

text: Thanks to “Rod Lm” I got some very nice input for the following problem. My intention is to simulate within a `Manipulate` function a given number of possible stock price paths and to calculate their mean. So, if something is manipulated the “new mean” is also calculated automatically and both are shown by a line and numerically within the chart. I would also like to know how I can also calculate the mean outside the chart. I tried the following but it didn't work:               Mean[Table[        RandomFunction[          GeometricBrownianMotionProcess[0.1, 0.2, 100], {0, 250, 0.05}][         "Path"], {20}][[All, -1]]]      The next step -if possible- is to introduce a lower boundary. So for example if the initial stock price is 100 then the boundary could be at, let’s say, 70. First, I want the boundary to be shown by a line. Second, I would like to stop all stochastic processes that fall below this threshold and then receiving the mean from all remaining processes. Alternatively, if it’s not possible to stop the processes that fall below the threshold, I also just want the mean of the remaining processes that are beyond the threshold. Third, I also would like to have a manipulable threshold. So, in short, I want to be able to switch the threshold between 20 and 90 and then receive the mean of the processes that remain beyond this boundary. Is something like that possible? If yes, I would be really thankful if I would receive some suggestions how to solve my problem, or at least parts of it. So far I have done the following and is it even possible to modify this code in such a manner that my problem can be solved? My code so far:                Manipulate[        SeedRandom[seed];         ListLogPlot[         Table[RandomFunction[          GeometricBrownianMotionProcess[μ,σ,S0], {0, 250, 0.05}]["Path"], {P}],          Joined -> True,          AxesLabel -> {"Time", "St"},          PlotLabel -> Style["Forecasted Stock Price\n (Brownian Motion)", Bold],          PlotRange -> All,          ImageSize -> 500,          PlotStyle -> Directive[{Thin, Lighter@Gray}]       ],        {{S0, 100, "Initial Stock Value"}, 100, 500, 0.05, Appearance -> "Labeled"},        {{μ, 0.01, "Drift μ"}, 0.01, 1, 0.05, Appearance -> "Labeled"},        {{σ, 0.01, "Standard Deviation σ"}, 0.01, 1, 0.05, Appearance -> "Labeled"},        {{P, 1, "Paths"}, 1, 100, 1, Appearance -> "Labeled"},        {{seed, 77777, "New Random Case"}, 10000, 999999, 1},       Button["Set Initial Values", {S0 = 25, μ = 0.01, σ = 0.01}, ImageSize -> 150],       ControlPlacement -> Left]      ![enter image description here](http://i.stack.imgur.com/jt1WL.png)


title: Solving differential equations with constraints

text: The problem is $$ x''(t)=F(t) - kx'(t) \\\ x'(0)=x(0)=0$$ with the additional constraints $$ 0 < x(t) <a \\\ 0 < x'(t) <b $$ where $F(t)$ is a known function and $a,b,k$ are given. At the moment I solve such problems with the help of _`WhenEvent[...]`_. **Example** :               First@First@NDSolve[{       x''[t] == 10 Sin[0.1 t] - 0.1 x'[t], x[0] == 0,        x'[0] == 0, WhenEvent[x[t] > 400, x[t] -> 399],        WhenEvent[x'[t] > 40, x'[t] -> 39],        WhenEvent[x[t] < 0, x[t] -> 1],        WhenEvent[x'[t] < -40, x'[t] -> -39]},       x, {t, 0, 200},       MaxSteps -> 50000];     Plot[{x[t] /. %, 10 x'[t] /. %}, {t, 0, 200}, PlotRange -> All,       PlotLegends -> {"x[t]", "x'[t]"}, ImageSize -> 700]      ![Plot of x, x'](http://i.stack.imgur.com/Mj9fD.png) What is the proper way to solve this type of problem?


title: Finding intersection of 3D lines

text: I want to find the intersection of 3D line. I want to detect the 3D light source position by intersecting 3D lines.               lightsource = {0, 0, 250};          p3d1 = {100, 100, 100};      p3d2 = {100, 0, 100};     p3d3 = {0, 100, 100};          d1={500/3, 500/3, 0};     d2={500/3, 0, 0};     d3={0, 500/3, 0};     (*equation of 3 ray*)     ray1 = Flatten[p3d1 + (d1 - p3d1)*t1];     ray2 = Flatten[p3d2 + (d2 - p3d2)*s1];     ray3 = Flatten[p3d3 + (d3 - p3d3)*k1];      I solved the system equation like this AX=B.               eq = ray1 - ray2 - ray3     var = Variables[eq]     A = Coefficient[#, var] & /@ eq;     B = eq - A.var     X = PseudoInverse[A].B (*in this case we can use Inverse[A] because A is invertible.*)     ray1 /. {t1 -> -1}     {100/3, 100/3, 200}(*It is not right*) **What is the problem here?**      ( _Second method: intersecting two rays by two rays_ )               eq = ray2 - ray3;     eqp = eq[[1 ;; 2]];     var = Variables[eqp];     coefs = N[Coefficient[#, var] & /@ eqp];     verif =       Chop[eq /. Solve[eqp[[1]] == 0 && eqp[[2]] == 0, var]]     sol = Solve[eqp == {0, 0}, var];     cs1 = Flatten[ray2 /. sol] (*It is true for the 3 rays*).          I am using this graphics                         Show[Graphics3D[Cuboid[{0, 0, 0}, {100, 100, 100}]],              Graphics3D[{Red, Thick, Line[{p3d1, d1}]}],              Graphics3D[{Green, Thick, Line[{p3d2, d2}]}],              Graphics3D[{Blue, Thick, Line[{p3d3, d3}]}]]      ![enter image description here](http://i.stack.imgur.com/admRs.png)


title: Writing an expression as sum of squares of expressions

text: Suppose we have a symmetric homogeneous polynomial expression $P$ in $X=(x_1,\cdots, x_n)$. I want to check whether there are functions $g(X)$ so that $P$ is of the form $\sum _{1\le i<j\le n} g_{ij}(X) (x_i-x_j)^2$. Is there a function to find $g$? As an example: $2 a^4 +2 b^4 + 2 c^4+ 4 a^3 b + 4 ab^3 + 4 a^3 c + 4 a c^3 + 4 b^3 c + 4 b c^3 - 4 a^2 b c- 4 ab^2 c - 4 a b c^2 \- 6 a^2 b^2 - 6 b^2 c^2 - 6 a^2 c^2= (a^2+b^2+2c^2+6ab)(a-b)^2+(b^2+c^2+2a^2+6bc)(b-c)^2+(c^2+a^2+2b^2+6ca)(c-a)^2 $ In _Mathematica_ :               2 a^4 + 2 b^4 + 2 c^4 + 4 a^3 b + 4 a b^3 + 4 a^3 c + 4 a c^3 +  4 b^3 c + 4 b c^3 -     4 a^2 b c - 4 a b^2 c - 4 a b c^2 - 6 a^2 b^2 - 6 b^2 c^2 - 6 a^2 c^2 ==     (a^2 + b^2 + 2 c^2 + 6 a b) (a - b)^2 + (b^2 + c^2 + 2 a^2 + 6 b c) (b - c)^2 +     (c^2 + a^2 + 2 b^2 + 6 c a) (c - a)^2


title: Simultaneously fitting multiple datasets

text: What is the proposed approach if one wants to simultaneously fit multiple functions to multiple datasets with shared parameters? As an example consider the following case: We have to measurements of Gaussian line profiles and we would like to fit a Gaussian to each of them but we expect them to be at the same line center, i.e. the fitting should use the same line center for both Gaussians. The solution I came up with looks a little clumsy. Any ideas on how to do this better, especially in cases where we have more than 2 datasets and more than one shared parameter? Example:               f[x_, amplitude_, centroid_, sigma_] :=        amplitude Exp[-((x - centroid)^2/sigma^2)]     data1 = Table[{x, RandomReal[{-.1, .1}] + f[x, 1, 1, 1]}, {x, -4, 6,               0.25}];     data2 = Table[{x, RandomReal[{-.1, .1}] + f[x, .5, 1, 2]}, {x, -8, 10,              0.5}];     gauss1 = NonlinearModelFit[data1, f[x, a1, x1, b1], {a1, x1, b1}, x,               MaxIterations -> 1000, Method -> NMinimize];     gauss2 = NonlinearModelFit[data2,               Evaluate[f[x, a2, x1, b2] /. gauss1["BestFitParameters"]], {a2,               b2}, x, MaxIterations -> 1000, Method -> NMinimize];     Join[gauss1["BestFitParameters"],gauss2["BestFitParameters"]]          datpl = ListPlot[{data1, data2}, Joined -> True,              PlotRange -> {{-10, 10}, All}, Frame -> True,              PlotStyle -> {Black, Red}, Axes -> False, InterpolationOrder -> 0];     Show[{datpl,            Plot[{Evaluate[f[x, a1, x1, b1] /. gauss1["BestFitParameters"]],                  Evaluate[                 f[x, a2, x1 /. gauss1["BestFitParameters"], b2] /.                  gauss2["BestFitParameters"]]}, {x, -10, 10}, PlotRange -> All,                  PlotStyle -> {Black, Red},                 Frame -> True, Axes -> False]}]      ![enter image description here](http://i.stack.imgur.com/V0H4f.jpg)


title: Fitting for multiple data in Mathematica

text: I've just started to use Mathematica for programming and I was wondering how I could do Fitting for several data sets. For example, I have one curve. For this curve I need to plot a fit curve, if my x>0, i'm using this equation: x a (1 - (0.2 - 0.5 x)/b), otherwise (if x<0) this one: x a (1 - (0.2 - 0.5 x)/b). a and b are two parameters I need to find to fit a curve. So to do this I used NonlinearModelFit: NonlinearModelFit[fl, If [x < 0.2, {x a (1 - (0.2 - 0.5 x)/b)}, {a (x - 0.04/(2 b))}], {a, b}, x]. And found my a and b. Is this correct? Now I have another curve and I have to do the same. But! This parameters (a and b) have to be the same for both curves. So how to plot two (three, four) fit curves for different data on one plot? Another problem is that for second curve an equation would be: If [x < 0.4, {x a (1 - (0.4 - 0.5 x)/b)}, {a (x - 0.16/(2 b))}], {a, b}, x] For third curve: If [x < 0.6, {x a (1 - (0.6 - 0.5 x)/b)}, {a (x - 0.36/(2 b))}], {a, b}, x]


title: Linking two differential equations

text: The solution and plot of the Differential equation is $\frac{d^2S(r)}{dr^2}+\frac{1}{r}\frac{dS(r)}{dr}-S(r)+3S^3(r)=0$               S0 = S[r] /.         NDSolve[{S''[r] + 1/r S'[r] - S[r] + S[r]^3 == 0,           S[1/10000000000000000] == 2.20620086,           S'[1/10000000000000000] == 0}, S, {r, 1/1000, 15}];     Plot[{S0}, {r, 1/100, 14}, PlotRange -> Full]      I have to use the above info in the differential equation: $\frac{d^2Z(r)}{dr^2}+\frac{1}{r}\frac{dZ(r)}{dr}-Z(r)+3S^2Z(r)-S^5(r)=0$ **And I have plotted like this: but it is not working**               sol = NDSolve[{s''[r] + 1/r s'[r] - s[r] + s[r]^3 == 0,          s[1/10000000000000000] == 2.20620086,          s'[1/10000000000000000] == 0}, s, {r, 1/1000, 15}];     sol = s[r] /. sol;     Z0 = Z[r] /.         NDSolve[{Z''[r] + 1/r*Z'[r] - Z[r] + 3 sol^2*Z[r] - sol^5 == 0,           Z[1/10000000000000000] == -16.17403,           Z'[1/10000000000000000] == 0}, Z, {r, 1/1000, 15}];     Plot[{Z0}, {r, 1/100, 14}, PlotRange -> Full]      The output curve will be like: ![enter image description here](http://i.stack.imgur.com/dFrMM.png). I have used r instead of $\rho$. I have tried several times to get the output but it shows some following which I don't understand. I have looked the error to get but couldn't. Sorry for the poor question. I have studied but couldn't fix it.


title: How to negate a boolean vector?

text: The `!True` returns `False` but how can I negate a boolean vector such as `{True, False, True}`? The exclamanation mark does not work here like this `!{True,False,True}`.


title: Newton-Raphson Method and the Van der Waal Equation Coding question

text: So I need to find Specific Volumes using the Newton Method from the Van der Waal equation over a change of constant Temperature but variant Pressure. I have been able to make a list of however many iterations of the altered Van der Waal equation for the root finding method from Pressure Min to Pressure Max (3.9atm to 59.2atm). Further, I am able to select one element/expression from that list and find the root/Specific Volume for that iteration. However, I cannot seem to figure out how to make a list of all the Specific Volumes, ie solve for the root in more than one expression at a time. This is my code solving for the 10th element in the list (about 12.2atm): * * * a = 1.36; b = .003183; R = .0820578; T = 333; inc = (59.2 - 3.9)/60; c = Table[ vi - ((((atm) + (a/vi^2)) (vi - b) - (R*T))/((atm) - (a/ vi^2) + (2 a*b/vi^3))), {atm, 3.9, 59.2, inc}] FindRoot[c[[10]] == 0, {vi, 4}, WorkingPrecision -> 20] * * * The output: * * *               {-((-27.3252 + (3.9 + 1.36/vi^2) (-0.003183 + vi))/(        3.9 + 0.00865776/vi^3 - 1.36/vi^2)) +        vi,           -((-27.3252 + (4.82167 + 1.36/vi^2) (-0.003183 + vi))/(        4.82167 + 0.00865776/vi^3 - 1.36/vi^2)) +        vi,           -((-27.3252 + (5.74333 + 1.36/vi^2) (-0.003183 + vi))/(        5.74333 + 0.00865776/vi^3 - 1.36/vi^2)) + vi,      ## AND SO ON............. Also: * * *               {vi -> 2.1930639624928306932}      * * * Which is a correct answer. I've tried this fix to solve for more than one element in the list(elements one to ten): * * *               FindRoot[c[[1;;10]] == 0, {vi, 4}, WorkingPrecision -> 20]      * * * but this yields an error: * * * > FindRoot::nveq: The number of equations does not match the number of > variables in FindRoot[c[[1;;10]]==0,{vi,4},WorkingPrecision->20]. >> * * * Any ideas? Thanks! EDIT: I'm sorry, it's working again for the Van der Waals but neither of the other 2. No matter how I change the starting point the output now is the same, odd number repeated and I don't know how to manipulate either FindMinimum nor Foldlist to change that. ## My code is               R = 8.314;     A = 507.28;     B = 0.10476;     a = 0.07132;     b = 0.07235;     c = 660000;     T = 333;     iter = 60;          FX = -P - (A (1 - a/v))/v^2 + (R T (1 - c/(T^3 v)) (B (1 - b/v) + v))/        v^2;     FXX = -((a A)/v^4) + (2 A (1 - a/v))/v^3 + (        R T (1 + (b B)/v^2) (1 - c/(T^3 v)))/v^2 + (        c R (B (1 - b/v) + v))/(T^2 v^4) - (        2 R T (1 - c/(T^3 v)) (B (1 - b/v) + v))/v^3;          FXXX = v - FX/FXX;          nn = Table[FXXX, {P, 400, 6000, iter}];          Ps = Table[P, {P, 400, 6000, iter}]          vs = Rest@FoldList[v /. Last@FindMinimum[#2, {v, #1}] &, 12, nn]          ListPlot[Thread[{Ps, vs}], Joined -> True, PlotRange -> {0, 7}]      * * * ## The output is               {1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77,       1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77,       1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77,       1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77,       1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77,       1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77,       1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77,       1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77,       1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77,       1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77,       1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77,       1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77,       1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77,       1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77,       1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77,       1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77,       1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77,       1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77,       1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77,       1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77,       1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77}      * * * I'm sorry I basically know nothing haha, but I need this done and I just don't know how to do it. Thanks so much for your help.


title: Boundary Condition Problem with Green's Function

text: My temperature distribution looks like:![enter image description here](http://i.stack.imgur.com/l06fu.png) Mean wire temperature is :               To[x_] := 2 c1 Cosh[x a1] + a2      And it is continuous at two points. In positive part:               Tu[x_] := 10 Exp[-m (x - 0.0005)]      In negative part, it is symmetry of Tu:               Tu2[x_] := 10 Exp[-n (- x - 0.0005)]      I found c1, m, n by using boundary conditions:               To[x_] := 2 c1 Cosh[x a1] + a2     Tu2[x_] := 10 Exp[-n (- x - 0.0005)]     bcs = {To[-L/2] == Tu2[-L/2], To'[-L/2] == a4 Tu2'[-L/2]};     s = First@Solve[bcs[[1]], c1]     Quiet@Solve[bcs[[2]] /. s[[1]], n]          To[x_] := 2 c1 Cosh[x a1] + a2     Tu[x_] := 10 Exp[-m (x - 0.0005)]     bcs = {To[L/2] == Tu[L/2], To'[L/2] == a4 Tu'[L/2]};     s = First@Solve[bcs[[1]], c1]     Quiet@Solve[bcs[[2]] /. s[[1]], m]      Now I know the constants of my functions. Then I have non-linear differential equation which can be solved with Green's function. I solved that differential equation with Green's function when the boundary conditions were To[L/2] ==0 and To[-L/2] ==0. Green's function must satisfy boundary condition and differential equation. Because of this reason               green[-(L/2)] == 0, green[L/2] == 0      Now, same problem but different boundary conditions which are presented above. Before, I used that code with To[L/2] ==0 and To[-L/2] ==0, and got green[y]               eqn1 = D[green[y], {y, 2}] q - green[y] (f0 - b g0 + I w p) - DiracDelta[x - y]     ans1 = DSolve[{eq1 == 0, green[-(L/2)] == 0, green[L/2] == 0}, green[y], y]      It worked well. But I could not applied new boundary conditions into those code. Program gave an error. Because I have 4 boundary conditions:               green[L/2] == Tu[L/2], green'[L/2] == a4 Tu'[L/2]     green[-L/2] == Tu2[-L/2], green'[-L/2] == a4 Tu2'[-L/2]      How can I get Green's function with new boundary conditions ?