0
$\begingroup$

I am new to Mathematica and am struggling with making plots of how features of solutions vary with parameters. For example, I can find the eigenvalues and zeroes of a Sturm Liouville problem (showing just the first nonzero eigenvalue for this example) with

sol = NDSolve[
{g''[x] - 
   4*c*(-((3 r sigma^3)/(L - x)^4) + (3 sigma^3)/(L + x)^4)*
    g'[x] + mu[x]*g[x] == 0, mu'[x] == 0,
 g[-b] == 1, g'[-b] == 0, g'[b] == 0},
{g, mu}, x,
Method -> {"EventLocator", "Event" -> g[x], 
  "EventAction" :> 
   Print["mu= ", lambda = mu[x], ",  g[", zero = x, "] = ", 
    g[x]]}
]

but I would like to see how these features vary with the parameters. Looking at the basic help files that come with Mathematica, I can loop over parameters using Do[(block of code given above),{c,cstart,cend,dc},{r,rstart,rend,dr}] etc. But then the documentation tells me how to make graphs by either typing thousands of numbers in curly braces or reading them from a file (of not obviously documented format).

I've spent two hours trying to store these numbers in arrays that ListPlot3d will recognize but I either get error messages when trying to create the arrays or error messages when ListPlot3d tries to read them. The closest I have come to success is writing the results to a sequential file

sfile = OpenWrite["zero.dat"]; Do[....
Write[sfile, {r, c, zero}]
....
Close[sfile]
z_list = Import["zero.dat", "Table"]

This gets the numbers in and out. If the imported data is displayed in table form it looks like a table of x,y,z data for ListPlot3d. But when I try to run it I get

In[55]:= ListPlot3D[z_list]
During evaluation of In[55]:= Rule::rhs: Pattern z_list appears on the right-hand side of rule System`ProtoPlotDump`modelData$26641[WrappedValues,z_list,Charting`Private`Tag$26644]->z_list. >>

During evaluation of In[55]:= ListPlot3D::arrayerr: z_list must be a valid array or a list of valid arrays. >>

Out[55]= ListPlot3D[z_list]

The idea that I have to write results to a file to plot them sounds like I should be invoking WATFOR. There must be an elegant way to plot the results of complex calculations, but in two hours looking through help (and google), I haven't found it. How do I do this?

$\endgroup$
  • $\begingroup$ What version of Mathematica do you have? Since version 9, there is ParametricNDSolve and WhenEvent has replaced "EventLocator". $\endgroup$ – Michael E2 Dec 21 '14 at 0:39
  • 1
    $\begingroup$ You shouldn't use _ in parameter names. What happens if you replace your z_list with, e.g., zList? $\endgroup$ – Karsten 7. Dec 21 '14 at 0:45
  • 1
    $\begingroup$ Or to interactively vary the parameters and see the effects, use Manipulate. reference.wolfram.com/language/ref/Manipulate.html $\endgroup$ – Bob Hanlon Dec 21 '14 at 0:54
  • $\begingroup$ Your code is confusing, but at the very least you can simplify it by replacing the constraint mu'[x] == 0 with mu == mu0 (a constant). Don't compound problems by worrying about how to read and write variables (Import and Export). Give us the simplest version of your code that works. $\endgroup$ – David G. Stork Dec 21 '14 at 1:30
  • 1
    $\begingroup$ Have you looked at the help file for NDSolve? There are many examples there that plot the output of the ODE. Maybe you could mimic one of these? $\endgroup$ – bill s Dec 21 '14 at 3:42
2
$\begingroup$

I woudl suggest defining your soulutions as a delayed function using SetDelayed (":=")

sol[c_, sigma_, L_, r_, b_] := ...

Now you can go ahead and generate a table of solutions like that:

sols = ParallelTable[sol[c, sigma, L, r, b], {c, cmin, cmax, dc}, {sigma, sigmamin, sigmamax, dsigma}, {L, Lmin, Lmax, dL}, {r, rmin, rmax, dr}, {b, bmin, bmax, db}]

and access a solution within sols by index such that

sols[[1,2,3,4,5]]

yields

sol[cmin + 0*dc, sigmamin + 1*dsigma, Lmin + 2*dL, rmin + 3*dr, bmin + 4*db]

you can also take subsets of the table you created and listplot them. I hope I understood your problem right.

PS: since your problem comes with 5 Parameters you get a 5 dimensional result wich is hard to plot, plotting subsets might be a solution. A good way to plot up to 3 dimensional results is using Image3D[].

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.