Community RSS Feed
http://community.wolfram.com
RSS Feed for Wolfram Community showing any discussions in tag Wolfram Scienceshowthread.php?s= sorted by activeSimplify Function
http://community.wolfram.com/groups/-/m/t/1475530
Hello,
I'm beginner with Mathematica and trying to get the basics.
I try to use the function called "Simplify".
I try to simplify a function Abs(Cos(x)) on 2 different intervals where the result is positive.
Thus, I would expect the Simplify function to provide Cos(x) as the result in both cases, exepect that it works inly for one of those intervals :
![enter image description here][1]
[1]: http://community.wolfram.com//c/portal/getImageAttachment?filename=SImplify2.jpg&userId=1475502
Do you please have an idea why ine the first case, result is not Cos[x] ?
Thanks for kind help
Best RegardsMorgan Bourgeois2018-09-23T14:23:59Zabout graphics error
http://community.wolfram.com/groups/-/m/t/1477208
I want to simulate Solar Planet's motion with mathematica, so I typed like this.
![The Code][1]
But there were some errors.
I can't figure out what is wrong and how to fix it.
![The Error][2]
[1]: http://community.wolfram.com//c/portal/getImageAttachment?filename=%EC%BA%A1%EC%B2%98.JPG&userId=1476693
[2]: http://community.wolfram.com//c/portal/getImageAttachment?filename=111.png&userId=1476693hojun lee2018-09-24T14:49:00ZThe MoleculeViewer package
http://community.wolfram.com/groups/-/m/t/1190717
I have just released a package entitled ["MoleculeViewer"][1], whose functionality is exactly what it says on the tin.
This package was inspired in part by previous efforts by [@BoB LeSuer][at0] and [@Bianca Eifert][at1]. I took the best parts of [their][2] [packages][3], along with some of the good parts of the built-in molecule renderer, and added a few of my own tweaks. One noticeable tweak would be the depiction of multiple bonds (just like what is done in some physical models), as in the following image:
MoleculeViewer["thiacloprid"]
![thiacloprid][4]
The package has a number of other nifty features and auxiliary functions, like highlighting:
MoleculeViewer["caffeine", Highlighted -> {"O", "N" -> Orange}]
![caffeine][5]
and legends:
MoleculeViewer[RunOpenBabel[GetChemSpider["calicheamicin", "InChI"]], PlotLegends -> True]
![calicheamicin][6]
Before using the package, you will need to install [Open Babel][7] for some of its conversion functionality. Additionally, to use the [ChemSpider][8] search functionality, you will need to [register][9] to obtain an API key.
Download the paclet from GitHub and [install in the usual manner][10]. Alternatively, using the technique featured [here][11], evaluate
PacletInstall["MoleculeViewer", "Site" -> "http://raw.githubusercontent.com/tpfto/MoleculeViewer/master"]
Documentation and a gallery are given as separate notebooks.
[at0]: http://community.wolfram.com/web/bobthechemist
[at1]: http://community.wolfram.com/web/biancaeifert
[1]: https://github.com/tpfto/MoleculeViewer/releases
[2]: https://github.com/biancaeifert/multi-bond-plot
[3]: https://github.com/bobthechemist/molviewer
[4]: http://community.wolfram.com//c/portal/getImageAttachment?filename=zotu8.png&userId=520181
[5]: http://community.wolfram.com//c/portal/getImageAttachment?filename=eTaCP.png&userId=520181
[6]: http://community.wolfram.com//c/portal/getImageAttachment?filename=nlZ63.png&userId=520181
[7]: http://openbabel.org/
[8]: http://www.chemspider.com/
[9]: http://www.chemspider.com/Register.aspx
[10]: https://mathematica.stackexchange.com/a/141888
[11]: https://mathematica.stackexchange.com/questions/155123J. M.2017-09-23T15:31:37Zslow program
http://community.wolfram.com/groups/-/m/t/1475465
Hi guys. I wrote a program to draw a plot. My program works for l=1 quickly. but for l>1 ( 2,3,...) it's so slow. I even wait a hour for ( l=3) but after that time it was still working! :( Please help me what can I do?reza hamzeh2018-09-23T15:58:19ZPJLink: Hooking up Mathematica and Python
http://community.wolfram.com/groups/-/m/t/1468475
Here's a cross-post of something I originally wrote [here](https://www.wolframcloud.com/objects/b3m2a1/home/pjlink-hooking-up-mathematica-and-python.html) about how to get python and Mathematica to work together like JLink.
I thought this might have broad appeal here.
PJLink: Hooking up Mathematica and Python
---
Mathematica is an incredibly powerful platform with a fun and intellectually pleasing language, but is incredibly expensive and closed source. Python is a convenient, pretty powerful language with a lot of support from the developer community. For as long as the two have existed people have been trying to tie them together, but very little has been done to do so at the native level with efficient, convenient exchange between the two. That's why over the past few weeks I took the time to build a clean, convenient link between the two. This post will go into how the link was built and some of its features, but first I think a little demo is appropriate.
## A Quick Demo
### Installing PJLink
The link is based off of the [J/Link](https://reference.wolfram.com/language/JLink/tutorial/Overview.html) interface built into Mathematica for hooking up Java and Mathematica. To wit, I called it [PJ/Link](https://github.com/b3m2a1/PJLink) . It lives on my paclet server as well as GitHub, so we can easily install it from there:
PacletInstall["PJLink", "Site"->"http://www.wolframcloud.com/objects/b3m2a1.paclets/PacletServer"]
(*Out:*)
![hookingupmathematicaandpython-274752603667507597](https://www.wolframcloud.com/objects/b3m2a1/home/img/hookingupmathematicaandpython-274752603667507597.png)
### Loading PJLink in Jupyter
For this demo we'll need the path to this thing as well (note that the version might change in the future):
%["Location"]
(*Out:*)
"~/Library/Mathematica/Paclets/Repository/PJLink-1.0.0"
Now we'll leave Mathematica and open up a Jupyter notebook:
![hookingupmathematicaandpython-921972091567300718](https://www.wolframcloud.com/objects/b3m2a1/home/img/hookingupmathematicaandpython-921972091567300718.png)
Next we'll get that path available so we can actually make use of the package. Then we'll load things from the subsidiary ```SubprocessKernel``` package which is included in the paclet and makes use of PJLink:
<pre>import os, sys
pjlink_path = &quot;~/Library/Mathematica/Paclets/Repository/PJLink-1.0.0&quot; #this is whatever path was extracted before
sys.path.insert(0, os.path.expanduser (pjlink_path))
from SubprocessKernel import SubprocessKernel
from SubprocessKernel import MathematicaBlock, LinkEnvironment
## these are helpers I&apos; ll use in the demo</pre>
![hookingupmathematicaandpython-695494794636171070](https://www.wolframcloud.com/objects/b3m2a1/home/img/hookingupmathematicaandpython-695494794636171070.png)
### Bidirectional Communication
Once we have this we can start a subprocess kernel which will open a Mathematica front-end to interact with. We'll also start and evaluator Mathematica can use to call back into python.
You may see a long string of output from your C compiler as the setup.py file builds out the native library that PJLink uses. Don't worry, this should only happen once. If it fails, raise an issue on [GitHub](https://github.com/b3m2a1/PJLink/issues) so I can deal with it.
Once Mathematica has loaded, we'll use the ```MathematicaBlock``` context manager so we can write something that looks a lot like Mathematica code and use the ```MEval``` function we'll define to run the code. That code for all this looks like:
<pre>ker = SubprocessKernel()
def MEval (expr, wait = True, kernel = ker) :
&quot;&quot; &quot;MEval evaluates a Mathematica expression in the Mathematica kernel
&quot; &quot;&quot;
kernel.drain() # just to make sure things are clen
return kernel.evaluate (expr, wait = wait)
ker.start()
ker.start_evaluator()</pre>
After that we can simply call into Mathematica:
<pre>with MathematicaBlock():
res = MEval (Set (M.hi, &quot;Hello from python!&quot;))
res</pre>
![hookingupmathematicaandpython-4447232891793681130](https://www.wolframcloud.com/objects/b3m2a1/home/img/hookingupmathematicaandpython-4447232891793681130.png)
We can see string ```"Hello from python!"``` was set to the symbol ```hi``` on the Mathematica side and was returned back by ```MEval``` . Symbols that aren't in the ```"System`"``` context need to be prefaced by an ```M.``` as that's a special class that can resolve symbol names like that.
We can also get efficient data transfer of arrays from either side. Here we'll take some Mathematica data and get it back out on the python side. The first thing we need to do is go to the Mathematica notebook that opened and load the ```"PJLink`"``` context. Then we'll install the python runtime that the ```SubprocessKernel``` object configured. This looks like:
<<PJLink`
InstallPython[ LinkObject->SubprocessKernel`$PyEvaluateLink, ProcessObject->None];
Once it's installed, we'll use it directly via ```PyEvaluate``` :
With[{arr= RandomReal[{-1, 1}, {50, 50, 50}]},
PyEvaluate[dat=arr]
]
Calls into python are done in an environment held only by the link, so to access that we need to wrap the evaluator we started ( ```ker.evaluator``` ) in a ```LinkEnvironment``` context manager:
<pre>with LinkEnvironment(ker.evaluator):
res = dat.shape
res</pre>
![hookingupmathematicaandpython-8804306407173974153](https://www.wolframcloud.com/objects/b3m2a1/home/img/hookingupmathematicaandpython-8804306407173974153.png)
Arrays are held as NumPy arrays by default on the python side, although this may be disabled. If disabled, they're held as a data type called ```BufferedNDArray``` which holds the data as a single C-contiguous array and allows slicing and viewing into it (although no efficient math or manipulation of any sort).
Finally, to close out the demo, we'll plot something on the Mathematica side and watch it come back on the python side. The code for this should be pretty self-explanatory by this point, but there is one cute feature to note:
<pre>with MathematicaBlock():
res = MEval(
Rasterize(
Plot(Sin (M.x), List (M.x, 0, Times (2, Pi)),
ImageSize = [250, 250],
PlotLabel = &quot;sin(x) as plotted in Mathematica&quot;
)
)
)
res</pre>
Unfortunately it really does matter that we pass a ```List``` expression instead of a python list for the second argument to ```Plot``` as otherwise the system hangs for reasons that aren't totally clear. On the other hand, we can see how nice options passing is in the interface. We make use of the python ```**kwargs``` setup and that ```ImageSize= ...``` and ```PlotLabel= ...``` both get automatically converted into rules (albeit with a ```String``` key instead of a ```Symbol``` ). The ```Rasterize``` is, sadly, similarly necessary as there is currently no logic in the package to automatically convert ```Graphics``` expression into their rasterized forms.
![hookingupmathematicaandpython-6192533434254394386](https://www.wolframcloud.com/objects/b3m2a1/home/img/hookingupmathematicaandpython-6192533434254394386.png)
I think we'll close out the demo here, though, and move onto a description of how this works.
## PJLink Native Library
The heart of PJLink is the C library that connects a python runtime to MathLink. The source for this can be found [here](https://github.com/b3m2a1/PJLink/blob/master/PJLink/PJLinkNativeLibrary/src/PJLinkNativeLibrary.cpp) . This library, once compiled by the setup.py file packaged with it, implements the basic MathLink calls in a way that python can use them and attempts to do so with efficient memory usage and data transfer.
### Data Sharing in the Native Library
The heart of the native library is the set of ```PutArray``` and ```GetArray``` functions it implements. Beyond anything else, it is the fast transfer of large arrays of data that makes a C-level connection so appealing. The way we handle this on the python side is via the python [buffer protocol](https://docs.python.org/3/c-api/buffer.html) . We enforce the condition that all data sent and received on the python side must be handled by an object that can work with a C-contiguous buffer of data. By default this is done with [NumPy](http://www.numpy.org/) if it is installed, but if not there is a custom object called ```BufferedNDArray``` in the [HelperClasses](https://github.com/b3m2a1/PJLink/blob/master/PJLink/HelperClasses.py) package that deals with this.
### Threading in the Native Library
Python has something called the [Global Interpreter Lock (GIL)](https://docs.python.org/3/c-api/init.html#thread-state-and-the-global-interpreter-lock) which is a method for synchronizing python state. Unfortunately for us, the presence of the GIL means that standard C calls of the kind we'll be using will cause all threads to lock. To get around this, every call into the MathLink library in the native library is wrapped in the ```MLTHREADED``` macro which handles the releasing and reacquiring of the lock. This allows our threads to work once more. Any extensions to the library should keep this in mind.
## Class Structures
PJLink provides a glut of classes that handle the details communication, so we will quickly detail what the important ones do. More information is always available upon request.
### The *Link classes
PJLink is based off of JLink and so it makes use of the same kind of class structure that JLink does. This means that it has a ```MathLink``` class that provides a template for the kind of link we'll work with and a ```KernelLink``` class that works specifically with Mathematica kernels. In general, we will only really work with a subclass of a ```KernelLink``` called a ```WrappedKernelLink``` that implements the ```KernelLink``` interface by calling into a ```NativeLink``` which is the only class which actually touches the native library at all.
If one is controlling a Mathematica kernel from python, it will be handled by a ```WrappedKernelLink``` .
### Reader class
The ```Reader``` class handles the other half of the communication. It waits for calls from Mathematica and processes them via the ```KernelLink._handlePacket``` function. Most commonly these calls in turn call ```KernelLink.__callPython``` which builds a python call from the symbolic python packet that ```PyEvaluate``` sends over the link. A ```Reader``` does its best not to completely prevent its link from passing data *to* Mathematica, but in general it is best not to depend on this as the ```NativeLink``` interface allows only a single thread to access the library at once for reasons of safety and stability.
### MathLinkEnvironment and MathLinkException
The ```MathLinkEnvironment``` is a standalone class that handles all of the various flags and state that the links need. It centralizes all information about what a given token or flag from MathLink means and provides utility functions for working with this. ```MathLinkException``` is a subclass of the standard python ```Exception``` class that handles the MathLink-specific exceptions that are returned. It in turn calls into ```MathLinkEnvironment``` to learn what various exceptions mean.
### MPackage, MLSym, and MLExpr
The HelperClasses package provides a large number of (generally) smaller classes that serve to make code cleaner in its implementation. A big part of this is done by the ```MPackage``` , ```MLSym``` , and ```MLExpr``` classes, which allow for a way to create packets with a syntax that looks more like standard Mathematica code. ```MLSym``` and ```MLExpr``` are types that a ```KernelLink``` knows how to put onto a link and ```MPackage``` provides utilities and a custom ```__getattr__``` so that the packet code can look like Mathematica code.
### MathematicaBlock and LinkEnvironment
Both ```MathematicaBlock``` and ```LinkEnvironment``` are also in the HelperClasses. They both edit the current evaluation state as [context managers](https://docs.python.org/3/reference/datamodel.html#context-managers) so that explicit references to ```MPackage``` can be dropped and variables held by a given link can be easily accessed. Being context managers, they are both used via ```with``` statements and change the execution environment of the enclosing block.
### BufferedNDArray, ImageData, and SparseArrayData
These are all data classes that allow for more efficient and convenient data transfer. The ```ImageData``` and ```SparseArrayData``` classes hold data coming in from Mathematica as put using ```PJLink`Package`AddTypeHints``` . They have methods to efficiently transform to more standard formats like ```PIL.Image``` and ```scipy.sparse.csr_matrix``` . As more data types are handled by ```AddTypeHints``` it can be assumed that more classes like these will be written.
## Mathematica-side Package
That was all to do with the python side of things, which is where most of the development work had to go. On the other hand, the Mathematica side of the equation still requires some explanation. The package itself is really quite simple, so please feel free to [peruse the source](https://github.com/b3m2a1/PJLink/blob/master/Mathematica/PJLink.wl) .
### InstallPython
```InstallPython``` is the most basic function in the package. It either finds or is given a python version or executable, attempts to open it via ```StartProcess``` , then links to it via ```LinkCreate``` and the ```start_kernel.py``` script provided in the PJLink python package.
Notably, all it really requires is a ```LinkObject``` , so you can pass one directly via the ```LinkObject``` option. It will also by default try to make a python ```ProcessObject``` but you can pass that via the ```ProcessObject``` option or you can pass ```None``` in which case it won't attach to a Mathematica controlled process.
### ClosePython
```ClosePython``` is the counterpart to ```InstallPython``` . It closes an opened python runtime by version or executable. When a new kernel is installed it is added to ```PJLink`Package`$PythonKernels``` and this is what ```ClosePython``` looks for to close.
### PyEvaluate / PyEvaluateString
This is the heart of the package. It takes Mathematica-esque code, converts it into a structure that can be processed by ```KernelLink.__callPython()``` and sends it over and waits for a response. The conversion is handled by ```PJLink`SymbolicPython`ToSymbolicPython``` which was originally written for the [PyTools package](https://github.com/b3m2a1/mathematica-PyTools) . This is the best way to move data to python as things like ```Image``` objects, packable arrays, and ```SparseArray``` objects will be moved over intelligently.
```PyEvaluateString``` is like ```PyEvaluate``` , but with the recognition that ```ToSymbolicPython``` will always be a little bit lacking. It allows one to simply call a string of python code on the link and get the results back.
### PyWrite / PyWriteString / PyRead / PyReadErr
These are all functions that make use of the fact that when the ```Reader``` object started it allowed an interactive session to keep running and reading / writing on stdin, stdout, and stderr. The ```Read``` functions read from stdout and stderr and the write functions write to stdin. The former takes Mathematica code and auto-converts it into a string. The latter simply passes in the given string.
## Future Work
PJLink 1.0.0, beefy as it already is, should only really be seen as the beginning. My hope is that much more can be done to allow for more native data type transfer and for intelligent communication between the two systems.
In my demo I tried to show some of the things that make the interoperation of the two so nice, but I obviously don't have the breadth of knowledge to know all of the many applications this can be put to. Applications built off of PJLink are always welcome and I'm happy to provide any requisite information and extensions to PJLink to get them built.
Alongside that, I think better integration on the Mathematica side is necessary. There is a partial interface for allowing a ```PythonObject``` structure to hide the details of ```PyEvaluate``` on the Mathematica side, but this needs work from both ends, first hooking up the ```Language`MutatationHandler``` interface and then extending the same on the python side. After that, a ```JavaBlock``` -like setup that allows a link to be opened, used, and cleaned up would be highly useful for sandboxing.
Finally, I'm sure there are numerous bugs hiding in the package as it stands. Please find them and let me know about them so they can be worked out.
In the meantime, I hope you enjoy PJLink and being able to use my two favorite languages symbiotically.b3m2a1 2018-09-20T09:18:51ZOptimizing release angle of a trebuchet (highschool project)
http://community.wolfram.com/groups/-/m/t/1474821
Hi, i'm new to using Mathematica and i'm trying to optimize the release angle of a trebuchet using Lagrangian multipliers. I have read much about it and I have come to understand how it is calculated. I ended up going with mathematica to solve the calculations since it is quite a task by hand. I've got a trebuchet with a counterweight pendulum, and a projectile on the end of a sling: ![enter image description here][1]
[1]: http://community.wolfram.com//c/portal/getImageAttachment?filename=Capture.PNG&userId=1474806
i am basing my calculations on the following info and equations https://classes.engineering.wustl.edu/2009/fall/ese251/presentations/(AAM_13)Trebuchet.pdf
I have attached my attempt in solving this using Mathematica. I have tested the method using a much simpler example with success(also attached). What I think is going wrong for me is the notation of the angular velocity and how i represent time.
I would like to find the angle(phi) and time where velocity(kinetic energy) of the projectile is the greatest with a potential energy constraint of 100.
My approach was to use
L(the,phi,t)=f(the,phi)−λ2(g1(the,phi)−c1)-λ2(g2(the,phi)−c2)
L(the,phi,t)=0
finding the gradient of
0=f(the,phi)−λ1(g1(the,phi)−c1)-λ2(g2(the,phi)−c2)
f(the,phi)=kinetic energy being optimized,
g1(the,phi)=potential energy,
c1=100,
g2(the,phi)=contreint for sling,
c1=0
Thanks in advancelasse thiellesen2018-09-23T13:03:35ZA Mathematica notebook that converts online lecture videos into LaTeX
http://community.wolfram.com/groups/-/m/t/1475783
Hi all, just sharing my project that was created at HackMIT 2018.
The premise of the project was to take online course videos (e.g. those on MIT OpenCourseWare), to feed it into a Mathematica notebook, which would then:
1. Find the key frames in the video where there is the most writing on the board
2. Break that frame down into lines of texts/equations in order
3. Send these cropped images to the API https://mathpix.com/ to convert them into LaTeX equations
Of course the resulting combined LaTeX would not be 100% accurate, but we believed this would save a lot of time in note-taking since modifying incorrect LaTeX is still much faster than LaTeX-ing from scratch. More details are available at https://devpost.com/software/r00t (note that you need your own API key for MathPix to function).
Some possible extensions that might be interesting to work on:
- Converting the Mathematica notebook to some kind of cloud service to be more easily accessible to the general public
- Allowing Mathematica to look up videos on YouTube automatically rather than needing to download the video prior to using this software
- A way to address professors' writing being slanted (up or down) rather than in a horizontal line: the current method of gaussian blurring in the horizontal direction to find individual equations would not work anymore
- Some form of machine learning to allow more accurate recognition of the equations written on the board
- Anything else you can think of that would be interesting to work on, really... :)
Thanks for reading!Jau Tung Chan2018-09-23T16:07:03ZWrong plot
http://community.wolfram.com/groups/-/m/t/1472475
Why it is showing plot like this?
![Wrong plot][1]
[1]: http://community.wolfram.com//c/portal/getImageAttachment?filename=ScreenShot2018-09-21at8.35.33PM.png&userId=1472442sivakkannan muthukumar2018-09-22T01:39:02ZPacking Multiple Ellipses in a Non-Convex Curve
http://community.wolfram.com/groups/-/m/t/1475425
![enter image description here][1]
[1]: http://community.wolfram.com//c/portal/getImageAttachment?filename=10ellipsesinoval.jpg&userId=29126
Introduction
Packing ellipses in a non-convex oval involves two problems:
keeping the ellipses inside the oval and keep the ellipses disjoint.
The basic approach for keeping the ellipses inside was discussed in a previous post http://community.wolfram.com/groups/-/m/t/1453823?p_p_auth=GJdd4EuK .
This function is applied to a group of discrete points,
from which an interpolation function is created for use in the optimization.
Keeping the ellipses disjoint, a simpler problem, is done by using Lagrange multipliers to create constraints, see http://www.optimization-online.org/DB_FILE/2016/03/5348.pdf .
The optimization is done as as random searches using ParametricIPOPTMinimize.
ParametricIPOPT Code
In[1]:= Needs["IPOPTLink`"]
In[2]:= statusrl = {0 -> "Solve_Succeeded", 1 -> "Solved_To_Acceptable_Level",
2 -> "Infeasible_Problem_Detected",
3 -> "Search_Direction_Becomes_Too_Small", 4 -> "Diverging_Iterates",
5 -> "User_Requested_Stop",
6 -> "Feasible_Point_Found", -1 -> "Maximum_Iterations_Exceeded", -2 ->
"Restoration_Failed", -3 -> "Error_In_Step_Computation", -4 ->
"Maximum_CpuTime_Exceeded", -10 -> "Not_Enough_Degrees_Of_Freedom", -11 ->
"Invalid_Problem_Definition", -12 -> "Invalid_Option", -13 ->
"Invalid_Number_Detected", -100 -> "Unrecoverable_Exception", -101 ->
"NonIpopt_Exception_Thrown", -102 -> "Insufficient_Memory", -199 ->
"Internal_Error"};
In[3]:=
In[4]:= iMin[obj_, cons_List, vwb_List, nrands_Integer, seed_Integer, opts___] :=
Block[{lecons, vars, lbubs, cfuns, clbubs, v, lb, ub, params, param, pip,
starts, pres, listres, status, goodres, finalres},
lecons = cons /. {
( a_ == b_ ) -> LessEqual[0, a - b, 0],
( a_ <= b_ ) -> LessEqual[-\[Infinity], a - b, 0],
(a_ >= b_ ) -> LessEqual[0, a - b, \[Infinity]]
};
{vars, lbubs} = Transpose[vwb /.
{v_, lb_?NumericQ, ub_?NumericQ} :> {v, {lb, ub}}];
If[Length[cons] == 0, {cfuns, clbubs} = { {}, {} },
{cfuns, clbubs} =
Transpose[lecons /. LessEqual[lb_, v_, ub_] -> {v, {lb, ub}} ]
];
params = Table[param[i], {i, Length[vwb]}];
pip = ParametricIPOPTMinimize[obj, vars, params, lbubs, cfuns, clbubs,
params, "RuntimeOptions" -> {"WarningMessages" -> False}, opts];
SeedRandom[seed];
starts = Transpose[RandomReal[#, nrands] & /@ lbubs];
pres = pip @@@ starts;
listres = {IPOPTMinValue[#], IPOPTArgMin[#], IPOPTReturnCode[#]} & /@ pres;
IPOPTDataDelete /@ pres;
status = listres[[All, 3]];
Print[Tally[status] /. {s_, n_Integer} :> { s /. statusrl, n}];
goodres = Select[listres, Or[#[[3]] == 0, #[[3]] == 1] &];
finalres = If[Length[goodres] == 0,
listres[[Ordering[listres[[All, 1]]]]][[1]],
goodres[[ Ordering[goodres[[All, 1]]]]][[1]]
];
{finalres[[1]],
Thread[vars -> finalres[[2]]], finalres[[3]] /. statusrl}
]
Packing Code
equation of circumscribing oval is oval[{x,y}] == 0
In[5]:= oval[{x_, y_}] = ((x - 1)^2 + y^2) ((x + 1)^2 + y^2) - (21/20)^4;
code to generate equation of an ellipse with semi-major axes "a" and "b"
In[6]:= rottransrl[{xc_, yc_, \[Theta]_}, {x_, y_}] =
Thread[{x, y} -> RotationMatrix[-\[Theta]].{x - xc, y - yc}];
In[7]:= Thread[{xel, yel} =
DiagonalMatrix[{a, b}^-1].RotationMatrix[-\[Theta]].{x - xc, y - yc}];
In[8]:= ell[a_, b_][xc_, yc_, \[Theta]_][x_, y_] = xel^2 + yel^2 - 1;
eliminate \[Lambda] from Lagrange multiplier equation for minimizing size of ellipse with axes "a" and "a/2" that overlaps oval
In[9]:= eld = Eliminate[
D[oval[{x, y}] == \[Lambda] ell[a, a/2][xc, yc, \[Theta]][x, y], {{x,
y}}], \[Lambda]] // FullSimplify
Out[9]= 5 y (xc + x (-2 + x xc) + xc y^2) +
3 (-xc y + (2 x - xc) y (x^2 + y^2) - x (-1 + x^2 + y^2) yc) Cos[
2 \[Theta]] ==
5 x (-1 + x^2 + y^2) yc +
3 (x^4 - x^3 xc + x (xc - xc y^2) - y (1 + y^2) (y - yc) +
x^2 (-1 + y yc)) Sin[2 \[Theta]] && a != 0
define function using previous result without a!=0
In[10]:= el[xc_, yc_, \[Theta]_] = eld[[1]]
Out[10]= 5 y (xc + x (-2 + x xc) + xc y^2) +
3 (-xc y + (2 x - xc) y (x^2 + y^2) - x (-1 + x^2 + y^2) yc) Cos[
2 \[Theta]] ==
5 x (-1 + x^2 + y^2) yc +
3 (x^4 - x^3 xc + x (xc - xc y^2) - y (1 + y^2) (y - yc) +
x^2 (-1 + y yc)) Sin[2 \[Theta]]
define function to find larger axis of largest ellipse centered at {xc,yc} with orientation \[Theta] that finds inside oval.
Min is needed due to multiple solution of Lagrange multiplier equations.
In[11]:= dist[{xc_?NumericQ,
yc_?NumericQ, \[Theta]_?NumericQ}] := -N @ Sign[oval[{xc, yc}]]*
Min[a /. NSolve[{a >= 0,
oval[{x, y}] == 0,
ell[a, a/2][xc, yc, \[Theta]][x, y] == 0,
el[xc, yc, \[Theta]]},
{a, x, y}, Reals]]
evaluate the function at certain discrete values
In[12]:= AbsoluteTiming[
Quiet[disttblflat =
N @ Flatten[
Table[{xc, yc, \[Theta], dist[{xc, yc, \[Theta]}]}, {xc, -2, 2, 1/
10}, {yc, -1, 1, 1/10}, {\[Theta], 0, 2 \[Pi], \[Pi]/10}], 2]];]
Out[12]= {797.969, Null}
construct an interpolation function from the values at discrete points
In[13]:= Dimensions[
disttbl = disttblflat /. {a_?NumericQ, b_, c_, d_} :> {{a, b, c}, d}]
Out[13]= {18081, 2}
In[14]:= itbl = Interpolation[disttbl, InterpolationOrder -> 2];
area of the oval for calculating packing fraction
In[15]:= ovalarea = NIntegrate[Boole[oval[{x, y}] <= 0], {x, -2, 2}, {y, -2, 2}]
Out[15]= 2.56423
plot of the oval
In[16]:= ovalplot =
ContourPlot[Evaluate[oval[{x, y}] == 0], {x, -1.5, 1.5}, {y, -1.5, 1.5},
ImageSize -> Small, ContourStyle -> Red];
forall[vars, c == d, a <= b] using Lagrange multipliers
In[17]:= lagforall[LessEqual[a_, b_], Equal[c_, d_], vars_List] :=
{Flatten @ {LessEqual[a, b], Equal[c, d], \[Lambda] >= 0,
D[a - b == \[Lambda] (c - d), {vars}]} /. \[Lambda] -> \[Lambda][
Sequence @@ vars],
Join[vars, {\[Lambda][Sequence @@ vars]}]}
constraint for two curves to not overlap
In[18]:= curvesnooverlap[curve1_, {xc1_, yc1_, \[Theta]1_},
curve2_, {xc2_, yc2_, \[Theta]2_}, {x_, y_}] :=
lagforall[
0 <= curve1[xc1, yc1, \[Theta]1][x, y],
curve2[xc2, yc2, \[Theta]2][x, y] == 0,
{x, y}
]
no overlap constraint for two ellipses
In[19]:= noeqs[i_, j_] =
curvesnooverlap[ell[a, a/2], {xc[i], yc[i], \[Theta][i]},
ell[a, a/2], {xc[j], yc[j], \[Theta][j]}, {xno[i, j], yno[i, j]}];
all the constraints: all the ellipses fit into the oval and all are disjoint
In[20]:= cons[n_Integer] :=
Join[Table[a == itbl[xc[i], yc[i], \[Theta][i]], {i, n}],
Flatten[Table[Sequence @@ noeqs[i, j][[1]], {i, n - 1}, {j, i + 1, n}], 1]]
all the variables with bounds
In[21]:= vwb[n_Integer] :=
Join[{{a, 0, 1}}, Table[{xc[i], -1.5, 1.5}, {i, n}],
Table[{yc[i], -.5, .5}, {i, n}],
Table[{\[Theta][i], \[Pi]/2, 3 \[Pi]/2}, {i, n}],
Flatten[Table[{xno[i, j], -1.5, 1.5}, {i, n - 1}, {j, i + 1, n}], 1],
Flatten[Table[{yno[i, j], -.5, .5}, {i, n - 1}, {j, i + 1, n}], 1],
Flatten[Table[{\[Lambda][xno[i, j], yno[i, j]], 0, 5}, {i, n - 1}, {j,
i + 1, n}], 1]]
plotting function, putting packing fraction over plot
In[22]:= plt[n_Integer, sln_] :=
Show[ovalplot,
ContourPlot[
Evaluate[Table[ell[a, a/2][xc[i], yc[i], \[Theta][i]][x, y] == 0, {i, n}] /.
sln[[2]]], {x, -1.5, 1.5}, {y, -1.5, 1.5}],
PlotLabel -> ToString[ ((n*\[Pi] a a/2)/ovalarea /. sln[[2]])]
]
Packing Calculation
In[23]:= plts = Range[10];
In[24]:= AbsoluteTiming[plts[[1]] = plt[1, iMin[-a, cons[1], vwb[1], 100, 0]]];
During evaluation of In[24]:= {{Maximum_Iterations_Exceeded,3},{Solve_Succeeded,97}}
In[25]:= AbsoluteTiming[
Do[plts[[i]] = plt[i, iMin[-a, cons[i], vwb[i], 50, 0]], {i, 2, 10}]]
During evaluation of In[25]:= {{Solve_Succeeded,40},{Maximum_Iterations_Exceeded,5},{Infeasible_Problem_Detected,5}}
During evaluation of In[25]:= {{Maximum_Iterations_Exceeded,26},{Infeasible_Problem_Detected,13},{Restoration_Failed,1},{Solve_Succeeded,10}}
During evaluation of In[25]:= {{Maximum_Iterations_Exceeded,31},{Infeasible_Problem_Detected,12},{Solve_Succeeded,7}}
During evaluation of In[25]:= {{Maximum_Iterations_Exceeded,36},{Solve_Succeeded,9},{Infeasible_Problem_Detected,5}}
During evaluation of In[25]:= {{Solve_Succeeded,11},{Maximum_Iterations_Exceeded,34},{Infeasible_Problem_Detected,5}}
During evaluation of In[25]:= {{Maximum_Iterations_Exceeded,35},{Solve_Succeeded,9},{Infeasible_Problem_Detected,6}}
During evaluation of In[25]:= {{Maximum_Iterations_Exceeded,43},{Solve_Succeeded,6},{Infeasible_Problem_Detected,1}}
During evaluation of In[25]:= {{Maximum_Iterations_Exceeded,35},{Solve_Succeeded,8},{Infeasible_Problem_Detected,7}}
During evaluation of In[25]:= {{Maximum_Iterations_Exceeded,36},{Solve_Succeeded,7},{Infeasible_Problem_Detected,7}}
Out[25]= {3989.62, Null}
In[26]:= pltsFrank Kampas2018-09-23T14:34:16ZTrying to create a table of dates
http://community.wolfram.com/groups/-/m/t/1472605
I was trying to create a table of the past 30 autumnal equinoxes, however I just can't seem to do it. If I try any thing fancy, Alpha just returns nearest autumnal equinox and that's it.
I've tried so far things like:
autumnal equinox [1988.. 2018.. 1]
and even went half-Mathematica with this one:
Table[autumnal equinox n, {n, 1988, 2018}]
Am I missing something? Am I taking too many liberties here? Please someone help.Poldek Novak2018-09-22T00:54:35ZHow we can delete same value from a Do loop?
http://community.wolfram.com/groups/-/m/t/1473026
Hi. I have a simple code. If I run the code the result are some numbers.
But, there are same numbers in the result. I want to remove those same numbers. Please help me. How can I do that?reza hamzeh2018-09-22T15:02:18ZWolfram Language / Mathematica Discord server chat (unoficial)
http://community.wolfram.com/groups/-/m/t/1472615
I have set up a Discord server for chatting with others about things related with the Wolfram Language / Mathematica, here: https://discord.gg/gj6VKAK
If you don't know about Discord, https://discordapp.com/Eric Parfitt2018-09-22T01:22:19ZWhy this Goto[end] statement does not work?
http://community.wolfram.com/groups/-/m/t/1472085
The Goto [end] statement does not function, and every time I run this program, it will show hold[Goto[end]]Jingxuan Sun2018-09-21T17:46:57ZOpenCVLink insight?
http://community.wolfram.com/groups/-/m/t/1472098
Recent versions of Mathematica contain OpenCVLink, presumably for support of image processing routines in WL. Its internal nature means that it is, understandably, relatively incomplete and undocumented.
I have some optic flow calculations to do, and, indeed, I see that two CUDA implementations are available,
$OpticalFlowFarneback
and
$OpticalFlowDualTVL1
but their parameters are not as easily decipherable as some of the other implementations (`$Dilation` for example).
I was wondering if anyone out there had successfully parsed the unparsable or, barring that, someone has a `.mc` link that uses OpenCV that I can just steal as a template.Maybe even WRI could think of releasing some of these incomplete undocumented links to the GitHub-o-verse where the public could enhance and refine them?Flip Phillips2018-09-21T18:08:51ZGet a numerical solution to a nonlinear ODE?
http://community.wolfram.com/groups/-/m/t/1458395
I am trying to solve a nonlinear ODE BY applying a NDsolve and using StiffnessSwitching method, but when I try to find the root of my equation it gives me a error message. The same code is working well in Mathematica version 9, but in M.version 11.3 that I just upgraded is not working I do not know why this happened.
Would anyone help me please?
Here is my code
Z=800;
g= 0.023800000000000000000;
k2= 0.000194519;
R= 1.5472;
ytest0= -13.911917694213733`;
ϵ = $MachineEpsilon ;
y1[ytest_?NumericQ] :=
NDSolve[
{y''[r] + 2 y'[r]/r == κ2 Sinh[y[r]] , y[1] == ytest,
y'[ϵ] == 0}, y, {r, ϵ, 1},
Method -> {"StiffnessSwitching", "NonstiffTest" -> False}];
y2[ytest_?NumericQ] :=
NDSolve[
{y''[r] + 2 y'[r]/r == κ2 Sinh[y[r]],
y[1] == ytest, y'[R] == 0}, y, {r, 1, R},
Method -> {"StiffnessSwitching", "NonstiffTest" -> False}];
y1Try[ytest_?NumericQ] := y'[1] /. y1[ytest];
y2Try[ytest_?NumericQ] := y'[1] /. y2[ytest];
f = ytest /. FindRoot[y1Try[ytest] - y2Try[ytest]==-Zg, {ytest, ytest0}]kolod al2018-09-15T03:44:54ZDerivative of Gaussian likelihood function?
http://community.wolfram.com/groups/-/m/t/1471082
Consider the following code:
In[63]:= f[x_]=1/(2*Pi*sigma^2)^0.5Exp[-(x-mu)^2/(2*sigma^2)]
In[64]:= p[x_]=Product[f[Subscript[x, i]],{i,1,n}]
In[68]:= D[p[x],mu]
![enter image description here][1]
It simply did nothing but added a partial derivative symbol at the front of p(x). Am I missing something here?
[1]: http://community.wolfram.com//c/portal/getImageAttachment?filename=Captured%E2%80%99e%CC%81cran2018-09-21a%CC%8003.35.21.png&userId=1471049Jason Wu2018-09-21T07:37:19ZEliminate G variable from this system of two equations?
http://community.wolfram.com/groups/-/m/t/1466372
Consider the following code:
Eliminate[{-A^16 + 3 A^11 F - 3 A^6 F^2 + A F^3 + 10 A^14 G -
20 A^9 F G + 10 A^4 F^2 G - 40 A^12 G^2 + 55 A^7 F G^2 -
15 A^2 F^2 G^2 + 80 A^10 G^3 - 85 A^5 F G^3 + 5 F^2 G^3 -
75 A^8 G^4 + 75 A^3 F G^4 - 25 A F G^5 + 75 A^4 G^6 -
75 A^2 G^7 + 25 G^8 + 25 A^6 H - 75 A^4 G H + 75 A^2 G^2 H -
25 G^3 H - 15 A^12 R + 5 A^7 F R + 10 A^2 F^2 R + 105 A^10 G R +
15 A^5 F G R + 5 F^2 G R - 300 A^8 G^2 R - 75 A^3 F G^2 R +
475 A^6 G^3 R + 25 A F G^3 R - 500 A^4 G^4 R + 375 A^2 G^5 R -
125 G^6 R - 25 A^8 R^2 + 25 A^3 F R^2 + 75 A^6 G R^2 +
50 A F G R^2 - 250 A^2 G^3 R^2 + 125 G^4 R^2 + 125 A^2 G R^3 ==
0 , -A^25 - 3125 A^10 B + 5 A^20 F - 10 A^15 F^2 + 10 A^10 F^3 -
5 A^5 F^4 + F^5 + 25 A^23 G + 15625 A^8 B G - 100 A^18 F G +
150 A^13 F^2 G - 100 A^8 F^3 G + 25 A^3 F^4 G - 275 A^21 G^2 -
31250 A^6 B G^2 + 850 A^16 F G^2 - 900 A^11 F^2 G^2 +
350 A^6 F^3 G^2 - 25 A F^4 G^2 + 1750 A^19 G^3 +
31250 A^4 B G^3 - 4000 A^14 F G^3 + 2750 A^9 F^2 G^3 -
500 A^4 F^3 G^3 - 7125 A^17 G^4 - 15625 A^2 B G^4 +
11375 A^12 F G^4 - 4500 A^7 F^2 G^4 + 250 A^2 F^3 G^4 +
19375 A^15 G^5 + 3125 B G^5 - 20000 A^10 F G^5 +
3750 A^5 F^2 G^5 - 35625 A^13 G^6 + 21250 A^8 F G^6 -
1250 A^3 F^2 G^6 + 43750 A^11 G^7 - 12500 A^6 F G^7 -
34375 A^9 G^8 + 3125 A^4 F G^8 + 15625 A^7 G^9 - 3125 A^5 G^10 +
25 A^21 R - 100 A^16 F R + 150 A^11 F^2 R - 100 A^6 F^3 R +
25 A F^4 R - 375 A^19 G R + 1125 A^14 F G R - 1125 A^9 F^2 G R +
375 A^4 F^3 G R + 2125 A^17 G^2 R - 4500 A^12 F G^2 R +
2625 A^7 F^2 G^2 R - 250 A^2 F^3 G^2 R - 4875 A^15 G^3 R +
6500 A^10 F G^3 R - 1500 A^5 F^2 G^3 R - 125 F^3 G^3 R -
1875 A^13 G^4 R + 3750 A^8 F G^4 R - 1875 A^3 F^2 G^4 R +
36250 A^11 G^5 R - 22500 A^6 F G^5 R + 1875 A F^2 G^5 R -
87500 A^9 G^6 R + 25000 A^4 F G^6 R + 103125 A^7 G^7 R -
9375 A^2 F G^7 R - 62500 A^5 G^8 R + 15625 A^3 G^9 R +
375 A^17 R^2 - 500 A^12 F R^2 - 125 A^7 F^2 R^2 +
250 A^2 F^3 R^2 - 6250 A^15 G R^2 + 6250 A^10 F G R^2 +
39375 A^13 G^2 R^2 - 25625 A^8 F G^2 R^2 + 1875 A^3 F^2 G^2 R^2 -
124375 A^11 G^3 R^2 + 48750 A^6 F G^3 R^2 - 2500 A F^2 G^3 R^2 +
215625 A^9 G^4 R^2 - 43750 A^4 F G^4 R^2 - 200000 A^7 G^5 R^2 +
12500 A^2 F G^5 R^2 + 75000 A^5 G^6 R^2 + 3125 F G^6 R^2 +
15625 A^3 G^7 R^2 - 15625 A G^8 R^2 - 1875 A^13 R^3 +
625 A^8 F R^3 + 1250 A^3 F^2 R^3 + 3125 A^11 G R^3 -
3125 A^6 F G R^3 + 25000 A^9 G^2 R^3 + 6250 A^4 F G^2 R^3 -
106250 A^7 G^3 R^3 - 3125 A^2 F G^3 R^3 + 175000 A^5 G^4 R^3 -
3125 F G^4 R^3 - 140625 A^3 G^5 R^3 + 46875 A G^6 R^3 -
3125 A^9 R^4 + 3125 A^4 F R^4 + 15625 A^7 G R^4 -
31250 A^5 G^2 R^4 + 31250 A^3 G^3 R^4 - 15625 A G^4 R^4 +
3125 A^5 R^5 == 0 }, G]mohamd fathi2018-09-19T11:28:30ZLazy lists in Mathematica
http://community.wolfram.com/groups/-/m/t/1467915
Hi all. In this post I want to demonstrate a package I wrote (and still tweak here and there), which implements Haskell-style lazy lists in Mathematica. Anyone who wants to play around with it can find here on Github:
https://github.com/ssmit1986/lazyLists
## What are lazy lists? ##
Before diving into implementation details, let's give some motivation for the use of lazy lists. A lazy list is a method to implicitly represent a long (possibly infinite) list. Of course, you cannot truly have an infinite list in your computer, so the central idea is to format the list as a linked list consisting of 2 parts. This means that a `lazyList` will always look like this:
lazyList[first, tail]
Here, `first` is the first element of the list and `tail` is a held expression that, when evaluated, will give you the rest of the list, which is again a `lazyList` with a first element and a tail. So in other words: elements of a `lazyList` are only generated when they are needed and not before. This makes it possible to represent infinite lists and perform list operations on them. To get the elements of the list, one can simply evaluate the tail as often as needed to progress through the list.
For example, let's define `lazyRange[]` as the lazy list of all positive integers. Then
Map[#^2&, lazyRange[]]
becomes the infinite list of all squares, again represented as a lazy list. You can go even further, though. For example, you can generate the triangular numbers by doing a `FoldList` over the integers and then select the odd ones with a `Select`:
Select[FoldList[Plus, 0, lazyRange[]], OddQ]
which is yet another lazy list. So if we want the first 100 odd triangular numbers, we simply evaluate the tail of this lazy list 99 times to get them. In contrast, if you'd try to do this with a normal list, you could do something like this:
Select[FoldList[Plus, 0, Range[n]], OddQ]
However, what value should you pick for `n`? If you pick it too low, you won't get your 100 numbers. If it's too high, you're doing too much work. Of course you could write some sort of `While` loop, but the code for that would be less concise and doesn't really play into the strengths of Wolfram Language.
## Implementation ##
To illustrate how my code works, I will reproduce some of the code in this blog post, though the package code is different in some respects for efficiency reasons.
The easiest way to prevent the tail from evaluating is to give `lazyList` the `HoldRest` attribute, which is how I implemented them:
Attributes[lazyList] = {HoldRest}
Next, we need some way to construct basic infinite lists like the positive integers. This is generally done recursively. My `lazyRange[]` function takes up to 2 arguments: a starting value (1 by default) and an increment value (also 1 by default):
lazyRange[start : _ : 1, step : _ : 1] := lazyList[start, lazyRange[start + step, step]]
We can extract the first element with `First` and advance through the list with `Last`:
First@lazyRange[]
First@Last@lazyRange[]
First@Last@Last@lazyRange[]
Out[100]= 1
Out[101]= 2
Out[102]= 3
We can also check that the tail of `lazyRange[]` is equal to the list of integers starting from 2:
In[103]:= Last@lazyRange[] === lazyRange[2]
Out[103]= True
Of course, iterating `Last` can be done with `NestList`, so if we want to get the first `n` elements of the lazy list, we can define the following special functionality for `Take` by setting an `UpValue` for `lazyList`:
lazyList /: Take[l_lazyList, n_Integer] := NestList[Last, l, n - 1][[All, 1]]
Take[lazyRange[], 10]
Out[105]= {1, 2, 3, 4, 5, 6, 7, 8, 9, 10}
As it turns out, nesting `Last` isn't actually the most efficient way to do this, so I ended up implementing `Take` with `ReplaceRepeated` and `Sow`/`Reap` to make the best use of the pattern matching capabilities of WL.
Next, we want to be able to do transformations on `lazyList`s. The simplest one is `Map`: you simply create a `lazyList` with the function applied to the first element and then `Map` the function over the tail:
lazyList /: Map[f_, lazyList[first_, tail_]] := lazyList[
f[first],
Map[f, tail]
];
Map[#^2 &, lazyRange[]]
Take[%, 10]
Out[117]= lazyList[1, (#1^2 &) /@ lazyRange[1 + 1, 1]]
Out[118]= {1, 4, 9, 16, 25, 36, 49, 64, 81, 100}
Similarly, `Select` is easily implemented by repeatedly evaluating the tail until we find an element that satisfied the selector function `f`. At that point we found our element and return a `lazyList`:
lazyList /: Select[lazyList[first_, tail_], f_] /; f[first] := lazyList[first, Select[tail, f]];
lazyList /: Select[lazyList[first_, tail_], f_] := Select[tail, f];
As an example, we can now find the first 10 numbers that are co prime to 12 and the first 10 squares that are 1 more than a multiple of 3:
Take[Select[lazyRange[], CoprimeQ[#, 12] &], 10]
Take[Select[Map[#^2 &, lazyRange[]], Mod[#, 3] === 1 &], 10]
Out[128]= {1, 5, 7, 11, 13, 17, 19, 23, 25, 29}
Out[129]= {1, 4, 16, 25, 49, 64, 100, 121, 169, 196}
I hope this gives a good enough overview of the benefits of `lazyLists` as well as giving you an idea of how to use them. In the package I tried to implement other list computational functionality (such as `MapIndexed`, `MapThread`, `FoldList`, `Transpose`, `Cases`, and `Pick`) for lazy lists as efficiently as possible.
Please let me know if you have further suggestions!Sjoerd Smit2018-09-19T22:04:27ZFunctions needed in Image Processing for future release of Mathematica
http://community.wolfram.com/groups/-/m/t/1467808
I am thrilled by the image processing capabilities built into Mathematica and witness its
power on a daily basis. I never have to use any other program for image processing,
save FIJI. In my opinion Mathematica's image processing core, albeit powerful and broad,
can adopt a thing or two from FIJI. As a biophysics student - who works with microscopy
images and other images in general - I firmly believe that the suggested functions mentioned
below can further add a lot of punch to the image processing core. Furthermore,
several of these functions have a role much broader than microscopy/fluorescence microscopy.
I sincerely hope that the image processing development team takes a note of the post.
Please let me know in the comments section if some function already exists. And feel free
to add functionalities to the list.
DESIRABLE FUNCTIONALITIES
- Hough Circle Transform (https://imagej.net/Hough_Circle_Transform)
- Creating surface through a 3D image by minimizing cost
(https://imagej.net/Minimum_Cost_Z_surface_Projection)
- Colocalization analysis
- Multiview reconstruction to form 3D images from different cameras (https://imagej.net/Multiview-Reconstruction)
- CLAHE (https://imagej.net/Enhance_Local_Contrast_(CLAHE))
- Rolling ball Background subtraction and/or better algorithms for background subtraction
- Scalebar and timestamp on images/videos (https://imagej.net/Time_Stamper)
- 16 bit image display in Mathematica
- PIV
- Kymograph
- Lookup Table to be applied to images
- images to FFT and FFT to images is much easier in FIJI (http://imagejdocu.tudor.lu/doku.php?id=gui:process:fft)
- Making ImageAlign robust to work on Binary masks and stack of images where objects are moving
- Bleach Correction for fluorescence microscopy (Mathematica's Histogram Transform can be extended)
(https://imagej.net/Bleach_Correction)
- Making superresolution images with CAIR (https://www.biorxiv.org/content/biorxiv/early/2017/12/19/236463.full.pdf)
- Introduction of an adjustable spline in manual mask creation tool (https://imagej.nih.gov/ij/docs/guide/146-27.html)
- Interactive affine transform/deformation of images
- Subpixel localization of small fluorescent particles in microscopy images
- Better tracking algorithms for cell and particle tracking for microscopy: my attempts in the direction are
(https://github.com/alihashmiii/Lineage-Mapper) and (https://github.com/alihashmiii/SMtrack . This can be made better
and sped up
[1]: https://imagej.net/Hough_Circle_TransformAli Hashmi2018-09-19T21:10:56ZInterpolate the following implicit function?
http://community.wolfram.com/groups/-/m/t/1466296
Hi,
I have a problem with solve. I need to define this function b of two variables, which is implicitly defined. I need it later to enter it into some other really nasty equations, which then again have to be solved numerically. So I wanted to define a grind of x and ss and then use an interpolation to get a smooth function version of b. However, I cannot even get a solution at a single point. (Maple does that without problems, so it shouldn't be a mathematical problem, but more of my stupidity). Below find the code an error text. If I use "NSolve" I even get no output. Any ideas?
Best
Christian
b[x_, ss_] := Solve[b == x - ss*Quantile[NormalDistribution, b], b, Reals]
b[0.3,0.1]
Solve::inex: Solve was unable to solve the system with inexact coefficients or the system obtained by direct rationalization of inexact numbers present in the system. Since many of the methods used by Solve require exact input, providing Solve with an exact version of the system may help.Christan Bauer2018-09-19T09:23:12ZIntegrate BesselJ function of zero order for galaxy?
http://community.wolfram.com/groups/-/m/t/1470004
Hi,
I trying to solve an integral to for a galaxy with parameters R (galaxy radius), r and z (cylindrical coordinates). The integral for f(r,z) is:
Integrate[Exp[-x z] BesselJ[0, R x] BesselJ[0, r x], {x, 0, Infinity}]
Mathematica can't solve this and I really don't know how to continue. Any suggestions on how to simplify the expression? Or how to solve it?
Thanks!Lina Persson2018-09-20T14:39:07ZHow to enforce priority on custom directory over paclets default repository
http://community.wolfram.com/groups/-/m/t/1471115
[**Topic cross posted in mathematica.stackexchange.com**](https://mathematica.stackexchange.com/q/182287/5478)
As explained in https://mathematica.stackexchange.com/q/59127/5478 loading packages from custom directories may not be so intuitive. Not to mention absolute lack of documentation of interaction between `$Path`, `PacletDirectoryAdd`, paclet's version etc.
Recent update to linked question explains something that hit me hard lately:
## The problem
If you have ``MyPackage` `` (V2) paclet installed in default repo and ``MyPackage` `` (V1) in customDir. Then:
PrependTo[ $Path, customDir];
PacletManager`PacletDirectoryAdd @ customDir;
<< MyPackage`
will load ``MyPackage` `` V2 from default repo.
I find this unexpected and really problematic
## Real world use case
You are working on dev branch on V2, part of testing involves packing it and installing locally, which means there will be V2 in default repo.
Before finishing V2 you need to push a fix to V1 so you create another fix branch from master (V1). You would like to load your code, **surprise** it loads V2 from default repo despite your changes to `$Path` and paclet directories.
## Questions
How to make sure that ``MyPackage` `` will be loaded from `customDir`?
``MyPackage` `` may contain 3rd part packages inside, e.g. `MyPackage/WL/GitLink` which it loads internally e.g. by temporarily blocking `$path` and prepending `MyPackage/WL`. However, `GitLink` may be already installed locally, newer or older, does not matter, the one from `customDir/MyPakcage/WL` should be loaded. Which means the solution can be based on fixing the issue for `MyPackage` only.
p.s. will try to provide a code to play with later.Kuba Podkalicki2018-09-21T07:00:42ZSolve analytically the following partial differential equations (PDE's)?
http://community.wolfram.com/groups/-/m/t/1395518
I have three PDE's and I want to solve it analytically. But I could not find any method to solve it. Can anyone suggest me which method is suitable for these type of PDE's.
Details are given in the attached file.
Anyone can help me it would be highly appreciated.Mirza Farrukh Baig2018-08-01T13:10:47ZChemical adsorption in fixed beds: PDE with boundary conditions
http://community.wolfram.com/groups/-/m/t/1470321
*NOTE: originally posted in response to:* http://community.wolfram.com/groups/-/m/t/1398247
----------
![Concentration Plot Zoomed][2]
# Related Notebooks
There are a couple of Mathematica notebooks related to chromatography that you may find useful. The first is from Professor Brian Higgins called [ChromatographyAnalysis.nb](https://sites.google.com/site/chemengwithmathematica/home/reaction-engineering/mathematica-notebooks-1/ChromatogaphyAnalysis.nb?attredirects=0&d=1) from his Ekaya Solutions website. The website has a plethora of Mathematica notebooks related to chemical engineering. If you have Wolfram SystemModeler, there is a [Chromatography Example Model](https://www.wolfram.com/system-modeler/examples/more/life-sciences/chromatography-column) from the Wolfram website.
# Basic Problem
To get the behavior that you desire, there is an implied step function that is missing from the equations that you posted. For "irreversible" adsorption, the time derivative of adsorbate loading must equal 0 at saturation. We also know that the solute concentration should be highest when the media is saturated. You will need to add a step function to drive the time derivative of W to be 0 at saturation. The result will be a phase boundary that moves and your equations are only valid for solid that is unsaturated. Therefore, the domain shape will be a function of time. I believe that this will mean that there are no analytical solutions available from DSolve and that you will need to use NDSolve to solve the system of PDEs. Traditionally, one would make a quasi-steady-state assumption and use a saturated moving boundary to obtain an approximate solution. I think an NDSolve solution will be instructive to show how the bed evolves over time and potentially give insights on simplified models.
# Dimensional Analysis
If you will indulge me for a moment, I would like to make the system of equations dimensionless so that we reduce the dependence to a minimum number of meaningful dimensionless groups with variables that have a nice scale.
There is a fluid phase and a solid phase in your system, so it makes sense that there should be a state variable for each phase. Our initial system of PDEs is given by:
$$\begin{gathered}
\varepsilon \frac{{\partial {c_d}}}{{\partial {t_d}}} + \left( {1 - \varepsilon } \right){\rho _p}\frac{{\partial {W_d}}}{{\partial {t_d}}} = - u\frac{{\partial {c_d}}}{{\partial {x_d}}} \qquad (1) \\
\left( {1 - \varepsilon } \right){\rho _p}\frac{{\partial {W_d}}}{{\partial {t_d}}} = {K_c}a\left( {{c_d} - c_d^*} \right) \qquad (2) \\
\end{gathered}$$
I added a subscript d to indicate that there are dimensions associated with those variables.
Next, to make variables dimensionless, simply divide the dimensioned variables by characteristic dimensions. For characteristic time, we will use the convective timescale or:
$$\begin{gathered}
{\tau _c} = \frac{{\varepsilon L}}{u} \\
t = \frac{{{t_d}}}{{{\tau _c}}};{t_d} = {\tau _c}t \\
x = \frac{{{x_d}}}{L};{x_d} = Lx \\
c = \frac{{{c_d}}}{{{c_0}}};{c_d} = {c_0}c \\
w = \frac{{{W_d}}}{{{W_{sat}}}};{W_d} = {W_{sat}}w \\
\end{gathered}$$
Now, let's return to equation (1) and make the appropriate substitutions.
$$\frac{\varepsilon }{{{\tau _c}}}{c_0}\frac{{\partial c}}{{\partial t}} + \left( {1 - \varepsilon } \right){\rho _p}\frac{{{W_{sat}}}}{{{\tau _c}}}\frac{{\partial W}}{{\partial t}} = - \frac{{u{c_0}}}{L}\frac{{\partial c}}{{\partial x}}\left\| {\frac{{\varepsilon L}}{u}} \right.\frac{1}{{\varepsilon {c_0}}}$$
Perform some algebra
$$\begin{gathered}
\frac{\varepsilon }{{{\tau _c}}}\frac{L}{{u{c_0}}}{c_0}\frac{{\partial c}}{{\partial t}} + \left( {1 - \varepsilon } \right){\rho _p}\frac{{{W_{sat}}}}{{{\tau _c}}}\frac{{\varepsilon L}}{u}\frac{1}{{\varepsilon {c_0}}}\frac{{\partial w}}{{\partial t}} = - \frac{{\partial c}}{{\partial x}} \\
\frac{{\partial c}}{{\partial t}} + \frac{{\left( {1 - \varepsilon } \right)}}{\varepsilon }\frac{{{\rho _p}}}{{{c_0}}}{W_{sat}}\frac{{\partial w}}{{\partial t}} = - \frac{{\partial c}}{{\partial x}} \\
m = \frac{{\left( {1 - \varepsilon } \right)}}{\varepsilon }\frac{{{\rho _p}}}{{{c_0}}}{W_{sat}} \\
\boxed{\frac{{\partial c}}{{\partial t}} + m\frac{{\partial w}}{{\partial t}} + \frac{{\partial c}}{{\partial x}} = 0} \\
\end{gathered}$$
The boxed formula represents the dimensionless version of (1). The dimensionless group, m, can be thought of as a loading factor (i.e., how much the adsorbant bed can adsorb versus the incoming concentration $c_0$. Now, we continue with equation (2).
$$\begin{gathered}
\left( {1 - \varepsilon } \right){\rho _p}\frac{{{W_{sat}}}}{{{\tau _c}}}\frac{{\partial w}}{{\partial t}} = {K_c}a{c_0}\left( {c - {c^*}} \right)\left\| {\frac{{{\tau _c}}}{{\varepsilon {c_0}}}} \right. \\
\frac{{\left( {1 - \varepsilon } \right)}}{\varepsilon }\frac{{{\rho _p}}}{{{c_0}}}{W_{sat}}\frac{{\partial w}}{{\partial t}} = \frac{{{K_c}a}}{\varepsilon }{\tau _c}\left( {c - {c^*}} \right) \\
m\frac{{\partial w}}{{\partial t}} = \frac{{{\tau _c}}}{{{\tau _a}}}\left( {c - {c^*}} \right);{\tau _a} = \frac{\varepsilon }{{{K_c}a}};n = \frac{{{\tau _c}}}{{{\tau _a}}} \\
\boxed{m\frac{{\partial w}}{{\partial t}} = n\left( {c - {c^*}} \right)} \\
\end{gathered} $$
The boxed formula represents the dimensionless version of (2). The new dimensionless group, n, can be thought of as the ratio of the convective time scale to the adsorption time scale.
Your system of PDEs for irreversible adsorption in dimensionless terms would be:
$$\begin{gathered}
\frac{{\partial c}}{{\partial t}} + m\frac{{\partial w}}{{\partial t}} + \frac{{\partial c}}{{\partial x}} = 0 \qquad (3) \\
m\frac{{\partial w}}{{\partial t}} = nc \qquad (4) \\
\end{gathered} $$
To correct equation (4) for the "equilibrium" condition of ${\left. {\frac{{\partial w}}{{\partial t}}} \right|_{w = 1}} = 0$ , we will need add a unit step function, $\theta$, to drive the time derivative to 0 at saturation and bringing over the left hand side of the equation or:
$$m\frac{{\partial w}}{{\partial t}} - nc\theta (1 - w)=0\qquad(4*)$$
# Bed Saturation Start-up Transient
If we start with a clean bed, it will take some time before the bed saturates at $x=0$. We can use equation (4) evaluated at $x=0$ to determine the time to initially saturate the front of the bed.
$$m{\left. {\frac{{\partial w}}{{\partial t}}} \right|_{x = 0}} = nc_0=n\qquad Conditions = \left\{ {\begin{array}{*{20}{c}}
{IC}&{w(t = 0,x = 0) = 0} \\
{FC}&{w(t = {t_{sat}},x = 0) = 1}
\end{array}} \right\}$$
$$w(t,x=0)=\frac{n}{m}t$$
$$\therefore t_{sat}=\frac{m}{n}$$
We can check our NDSolve solutions to verify that they approximate this relation.
# Mathematica Implementation
We should be able to use NDSolve to solve this system of PDEs. I had some performance issues with UnitStep\[\] so I opted to use the Tanh\[\] function to create a sharp, but differentially continuous approximation to the UnitStep\[\] function. Perhaps, there is a more elegant Mathematica implementation, but the Tanh\[\] functions seems to be working for $m$ and $n$ < 20.
Finally, we need two initial conditions (the bed and liquid is initially free of adsorbate species) and one boundary condition (a "UnitStep" of species at the inlet).
The initial conditions:
$$\begin{gathered}
c[0,x] = 0 \\
w[0,x] = 0 \\
\end{gathered} $$
The boundary condition:
$$c[t,0] = \tanh (200t)$$
This yields the following set of equations.
eqn1 = D[cp[t, x], t] + m D[wp[t, x], t] + D[cp[t, x], x] == 0;
eqn2 = m D[wp[t, x], t] -
n cp[t, x] (1 - (1 + Tanh[200 (wp[t, x] - 1)])/2) == 0;
ics = {wp[0, x] == 0, cp[0, x] == 0};
bcs = {cp[t, 0] == Tanh[200 t]};
eqns = {eqn1, eqn2} ~ Join ~ ics ~ Join ~ bcs;
opts = Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 5000}};
The dimensional analysis says that this system of equations depends on only two dimensionless parameters suggesting that ParametricNDSolveValue\[\] might be the flavor of NDSolve\[\] to use. To resolve the sharp features, we will specify a large number of spatial points as shown in opts and a small timestep as shown below to return a parameteric function, _pfun_.
pfun = ParametricNDSolveValue[
eqns, {cp, wp}, {t, 0, 10}, {x, 0, 1.1}, {m, n}, opts,
MaxStepFraction -> 1/1000];
With a parametric function, it is easy to return an interpolation function, _ifun_, based on the parameters $m$ and $n$. We will set the loading factor, $m$, to 5 and the convective timescale to adsorption timescale, $n$, to 10 (we want the adsorption to be fast relative to convection). The following takes a couple of minutes to run on my machine.
ifun = pfun[5, 10];
c = ifun[[1]];
w = ifun[[2]];
We can plot the particular solution of the fluid solute concentration, $c(t,x)$, with Plot3D\[\].
Plot3D[Evaluate[{c[t, x]}], {x, 0, 1}, {t, 0, 10},
PlotRange -> All, PlotPoints -> {200, 200},
ColorFunction -> (ColorData["DarkBands"][#3] &),
MeshFunctions -> {#2 &}, Mesh -> 20, AxesLabel -> Automatic,
MeshStyle -> {Black, Thick}, ImageSize -> Large]
![Concentration Plot][1]
We see some discontinuities at small $t$ as the bed starts loading to saturation. We can zoom into the small $t$ and $x$ for a clearer picture.
Plot3D[Evaluate[{c[t,x]}], {x, 0, 0.35}, {t, 0, 1},
PlotRange -> All, PlotPoints -> {200, 200},
ColorFunction -> (ColorData["DarkBands"][#3] &),
MeshFunctions -> {#2 &}, Mesh -> 50, AxesLabel -> Automatic,
MeshStyle -> {Black, Thick}]
![Concentration Plot Zoomed][2]
We predicted that the bed at $x=0$ should staturate at $t=\frac{m}{n}=\frac{5}{10}=0.5$ and the plot does indicate that the bed does indeed begin to saturate at the predicted time. We also see that the bed saturation front grows linearly in time. Also, post saturation, the concentration falls off approximately exponentially with $x$ at constant time from the saturation boundary. Prior to initial saturation, the MeshFunctions curves do not appear that they can be described by simple functions. Even if a DSolve solution existed, it probably would not have a simple form.
A countour plot of $w$ and $c$ after initial saturation may give us some insights to a simplified model. The first contour plot is of the solid loading.
ContourPlot[
Evaluate[{w[t + m/n, x]}], {x, 0, 5/(m + 1.3)}, {t, 0, 5},
ColorFunction -> "DarkBands", Mesh -> 40, MeshFunctions -> {#3 &},
FrameLabel -> {"x", "t"}, ImageSize -> Large]
![Solids Loading][3]
A remarkably linear _Saturation Boundary_ forms after $t=\frac{m}{n}$. The contour lines look essentially parallel implying that we might model $w$ as function of a single variable perpendicular to the phase boundary. The concentration countour plot looks similar but there is an accumulation error that grows linearly in time.
ContourPlot[
Evaluate[{c[t + m/n, x]}], {x, 0, 5/(m + 1.3)}, {t, 0, 5},
ColorFunction -> "DarkBands", Mesh -> 40, MeshFunctions -> {#3 &},
PlotRange -> {0, 1}, FrameLabel -> {"x", "t"}, ImageSize -> Large]
![Fluid Concentration Contour Plot][4]
If we restrict the domain to only consider the unsaturated bed, then we can eliminate $w$ by combining equations (3) and (4) because the UnitStep function containing $w$ can be ignored.
$$\begin{gathered}
\frac{{\partial c}}{{\partial t}} + m\frac{{\partial w}}{{\partial t}} + \frac{{\partial c}}{{\partial x}} = 0 \qquad (3) \\
m\frac{{\partial w}}{{\partial t}} = nc \qquad (4) \\
\\
\frac{{\partial c}}{{\partial t}} + \frac{{\partial c}}{{\partial x}}+ nc = 0 \qquad (5)
\end{gathered} $$
You can use Mathematica to evaluate the partial derivatives somewhere in the unsaturated region to assess their magnitudes.
D[c[t, x], t] /. {t -> 2.5, x -> 0.4} (*0.8176759782035153`*)
D[c[t, x], x] /. {t -> 2.5, x -> 0.4} (*-4.995343977849113`*)
Mathematica suggests that the spatial derivative is much larger in absolute magnitude than the time derivative. We could make the pseudo steady-state assumption and (5) becomes a simple first order ODE as shown in (6).
$$\frac{{\partial c}}{{\partial x}} + nc = 0 \qquad (6) $$
More rigorously, but still approximately, we could try a transformed coordinate system as depicted in the diagram below so that the moving boundary line given by $t=m x + \frac{m}{n}$ coincides with the ${t}'$ axis. In this reference frame, (5) should depend only on one independent variable ${x}'$ and the moving boundary goes away.
![Transformed Boundary][5]
To transform $(x,t)\rightarrow ({x}',{t}')$ we can use Mathematica's Transform functions to first rotate the original coordinate system by $-\theta=ArcTan[m,1]$ then translate up to $\frac{m}{n}$. Transformations can always be confusing so it is good to have a check to see if the transform is valid. In the transformed coodinates, the moving boundary line should simply be ${x}'=0$ .
theta = ArcTan[m, 1];
rt = RotationTransform[-theta];
tt = TranslationTransform[{0, m/n}];
trf = tt@*rt;
tr = trf@{xp, tp} // Simplify;
rules = Thread[{x, t} -> tr];(*{x -> (tp + m*xp)/Sqrt[1 + m^2], t -> m*(n^(-1) + tp/Sqrt[1 + m^2]) - xp/Sqrt[1 + m^2]}*)
We can apply the transformation rules to the moving boundary line to show that indeed ${x}'=0$.
(t == m x + m/n) /. rules // Simplify (*Sqrt[1 + m^2]*xp == 0*)
But we want the inverse transform to get ${x}'$ in terms of $x$ and $t$, which we can obtain with the InverseFunction\[\].
invtrf = InverseFunction[trf];
invtr = invtrf@{x, t} // Simplify;
invrules = Thread[{xp, tp} -> invtr] // InputForm
(*{xp -> (m - n*t + m*n*x)/(Sqrt[1 + m^2]*n), tp -> (-m^2 + m*n*t + n*x)/(Sqrt[1 + m^2]*n)}*)
Revisting (5), we will combine the $x$ and $t$ variables into one independent variable ${x}'$ using the chain rule as shown below:
\begin{gathered}
\frac{{\partial c}}{{\partial t}} + \frac{{\partial c}}{{\partial x}} + nc = 0; \\
\frac{{\partial x'}}{{\partial t}}\frac{{\partial c}}{{\partial x'}} + \frac{{\partial x'}}{{\partial x}}\frac{{\partial c}}{{\partial x'}} + nc = 0; \\
\left( {\frac{1}{n}} \right)\left( {\frac{{\partial x'}}{{\partial t}} + \frac{{\partial x'}}{{\partial x}}} \right)\frac{{\partial c}}{{\partial x'}} + c = 0 \\
\end{gathered}
Mathematica can evaluate and simplify the change of variables derivatives as shown below.
1/n (D[#, t] + D[#, x]) & [xp /. invrules] // Simplify (*(-1 + m)/(Sqrt[1 + m^2]*n)*)
Leaving us with equation (7) and the corresponding solution (we could use DSolve, but the solution to this ODE is straightforward).
$$\begin{gathered}
\frac{{m - 1}}{{n\sqrt {1 + {m^2}} }}\frac{{dc}}{{dx'}} + c = 0 \qquad (7) \\
c({x}')=e^{-\frac{n\sqrt{1+m^{2}}}{m-1}{x}'} \\
\end{gathered}$$
We can substitute and simplify $x$ and $t$ for ${x}'$ using the $invrules$ previously defined.
E^(-((Sqrt[1 + m^2]*n*xp)/(-1 + m))) /. invrules // InputForm (* E^(-((m - n*t + m*n*x)/(-1 + m))) *)
We recognize that $c$ is only defined to the right of the moving boundary to obtain the following piecewise function.
$$c(t,x)=\begin{Bmatrix}
e^{\frac{n \left(\left(t-\frac{m}{n}\right)- m x\right)}{m-1}} & m x>t-\frac{m}{n}\\
1 & \text{True}
\end{Bmatrix}$$
If you wanted to wrap the whole process in one Mathematica call, you could do something like the following.
cc[xp /. invrules] /. First@DSolve[{
(1/n (D[#, t] + D[#, x]) & [xp /. invrules]) D[cc[xp], xp] +
cc[xp] == 0,
cc[0] == 1},
cc,
xp
] // InputForm (* E^(-((m-n*t+m*n*x)/(-1+m))) *)
We will now create a function called _breakthrough_ for the simplified model and a plotting function to compare the simple model to the result obtained by ParametricNDSolveValue\[\].
(* Simplified model *)
breakthrough[m_, n_, t_, x_] :=
Piecewise[{{Exp[n ((t - m/n) - m x)/(m - 1)], m x > (t - m/n)}, {1,
True}}]
(* Plot function to compare simplified model to NDSolve for given m,n, and t *)
pltfn := Plot[{Callout[breakthrough[#1, #2, #3, x], "Model", (
Log[2] (-1 + #1) - #1 + #2 #3)/(#1 #2),
CalloutMarker -> "CirclePoint"],
Callout[c[#3, x], "NDSolve", (
Log[2] (-1 + #1) - #1 + #2 #3)/(#1 #2),
CalloutMarker -> "CirclePoint"]}, {x, 0, 1},
PlotRange -> {0, 1.05}, Filling -> {1 -> {2}}, Exclusions -> None,
PlotStyle -> Thick, Frame -> True,
FrameLabel -> {Style["x", 18], Style["c(t,x)", 18],
Style[StringForm["m=``, n=``, and t=``", #1, #2,
NumberForm[#3, {2, 1}, NumberPadding -> {" ", "0"}]], 18]},
ImageSize -> Large] & ;
(* Create a series of images at different times for animation *)
plts := Table[pltfn[#1, #2, t], {t, #1/#2, 10, 0.1}] &;
Here is an animation of the current solution for $m=5$ and $n=10$.
Export["510.gif", plts[5,10], "AnimationRepetitions" -> \[Infinity]]
![m=5 n=10 animated gif][6]
The accumulation term causes NDSolve to lag our simple model. At larger $m$ and $n$ values, this lag is reduced. It is easy to insert $m=20$ and $n=10$ into our parametric function and animate the result to verify the accumulation lag is reduced. Now, the simple model approximates well the NDSolve solution.
ifun = pfun[20, 10];
c = ifun[[1]];
Export["2010.gif", plts[20, 10],
"AnimationRepetitions" -> \[Infinity]]
![m=20 n=10 animated gif][7]
# Summary
- One needs to add a "UnitStep" function to stop adsorption when the bed becomes saturated.
- This will create a moving boundary meaning that the shape of the domain changes with time.
- Probably means DSolve will not find a solution.
- Dimensional analysis simplifies the equations and reduces the number of parameters that the equations depend on.
- ParametricNDSolveValue\[\] can be used to solve the system of PDEs.
- Using Plot3D and ContourPlot, we gain insights on a simplified model for the moving boundary.
- We can use Mathematica's Transform functions to create a new coordinate system that is perpendicular to the moving saturation boundary creating a simple solvable ODE.
- We created a simple model that tracks the results of NDSolve and the model should only get better at larger $m$ and $n$ values.
[1]: http://community.wolfram.com//c/portal/getImageAttachment?filename=ConcentrationPlot.png&userId=1402928
[2]: http://community.wolfram.com//c/portal/getImageAttachment?filename=ConcentrationPlotStartUp.png&userId=1402928
[3]: http://community.wolfram.com//c/portal/getImageAttachment?filename=LoadingContour.png&userId=1402928
[4]: http://community.wolfram.com//c/portal/getImageAttachment?filename=ConcentrationContour.png&userId=1402928
[5]: http://community.wolfram.com//c/portal/getImageAttachment?filename=SaturationBoundary.png&userId=1402928
[6]: http://community.wolfram.com//c/portal/getImageAttachment?filename=510.gif&userId=1402928
[7]: http://community.wolfram.com//c/portal/getImageAttachment?filename=2010.gif&userId=1402928Tim Laska2018-09-20T18:39:38ZSolve a PDE with boundary conditions (chemical adsorption in fixed beds)?
http://community.wolfram.com/groups/-/m/t/1398247
Dear Wolfram team:
I have been trying for week to solve a system of 2 partial differential equations describing the adsorption of a chemical substance on a fixed bed (for example, a column of activated carbon). The 2 equations are the following, taken from McCabe (1993):
![Description of eq 1][1]
![Description of eq 2][2]
[1]: http://community.wolfram.com//c/portal/getImageAttachment?filename=EQ1.png&userId=1020580
[2]: http://community.wolfram.com//c/portal/getImageAttachment?filename=EQ2.png&userId=1020580
Unfortunately I cannot get past the general solution (with arbitrary constants) because when I try to put boundary conditions the Mathematica program fails. Maybe I am using the wrong command or syntax, or maybe there are too much or too few boundary conditions.
I have left attached the program, where I tryed to simplify the problem combining both equations in a third.
Thank you in advance for your help.
Best regards,
Alberto SilvaAlberto Silva Ariano2018-08-06T01:04:00ZCoupled convective diffusive PDE in multiple regions with discontinuities
http://community.wolfram.com/groups/-/m/t/1470252
*NOTE: originally posted in response to:* http://community.wolfram.com/groups/-/m/t/1395518
----------
![enter image description here][10]
I looked at this problem as an opportunity to build some modeling and simulation skills within the _Mathematica_ framework and to brush up on an area of physics that I have not looked at in decades. Solving a coupled convective diffusive PDE in multiple regions with discontinuous physical properties can be quite challenging. In the end, I was pleasantly surprised that I was able to "model" this system and obtain reasonable results without too much difficulty. I am documenting my path here in the hopes that someone will find it useful.
# Preliminaries
I agree that there is no analytical solution to this system of coupled PDEs, and that a numerical solution is in order. The solution to the simplest convection-diffusion equation $$\rho C_pv_x\frac{\partial T}{\partial x}=k \frac{\partial^2 T}{\partial y^2}$$ for a parallel plate channel with fully developed uniform flow subjected to a uniform heat flux, [Sparrow and Siegel, 1960, pg 167](https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/20150018974.pdf) showed is an infinite series. For a coupled problem, I would presume that the infinite series coefficients would depend on the infinite series coefficients of the other variables making it intractable. [Ouyang _et al_](http://vafai.engr.ucr.edu/Publications/2013a/Ouyang.pdf) obtain an approximate solution with an infinite series solution to a thermally developing Darcy flow in porous media.
In the FEM, it is very easy to over constrain or under constrain a system leading to bizarre results. The developers of particular tools usually have developed a framework that facilitates the construction and solving of models. It is best to operate within that framework so that we can create a functioning model with least amount of effort.
Also, in the FEM, the mesh is critcal to finding an accurate solution. In your case, you have two regions with distinct physical properties. We expect that there will be extreme changes in behavior at the interface of the two regions so we will want to resolve that area with a finer mesh.
Finally, I would recommend to start small and build complexity verifying and validating along the way.
# System Description
I always like to begin with a picture of the system I am trying to model. My understanding is that there is a fluid flowing over a porous material as shown in the picture below. The flow in the x direction is assumed to be laminar parabolic flow in the fluid domain and constant uniform in the porous domain.
#### System Sketch
![enter image description here][1]
Physically, there must be an impermeable solid barrier between the two domains otherwise we would expect some flow in the Y direction. We will take advantage of this fact later.
# Porous Media Domain Only
### Porous Media Energy Transport Equations
As stated previously, the developers have created a framework to facilitate FEM modeling and simulation. For me, I find it often is easier to start over with that in mind than to start in the middle.
The equations below for energy transport with Local Themal Non-Equilibrium (LTNE) look formally similar to yours if we eliminate the transient and volumetric source terms shown in $\color{Red}{Red}$. If they are not, you should be able to follow the logic to develop your own.
$$\begin{gathered}
{\color{Red}{(1 - \varepsilon ){\left( {\rho {{\hat C}_p}} \right)_s}\frac{{\partial {T_s}}}{{\partial t}}}}+ {\left( {\rho {{\hat C}_p}} \right)_s}\left( {{\mathbf{v}}_s \cdot \nabla {T_s}} \right) - (1 - \varepsilon )\nabla \cdot {k_s}\nabla {T_s} - {\color{Red}{(1 - \varepsilon )q_s^{'''}}} - h({T_f} - {T_s}) = 0 \qquad (1*) \\
{\color{Red}{ \varepsilon {\left( {\rho {{\hat C}_p}} \right)_f}\frac{{\partial {T_f}}}{{\partial t}}}} + {\left( {\rho {{\hat C}_p}} \right)_f}\left( {{\mathbf{v}} \cdot \nabla {T_f}} \right) - \varepsilon \nabla \cdot {k_f}\nabla {T_f} - {\color{Red}{\varepsilon q_f^{'''}}} - h({T_s} - {T_f}) = 0 \qquad (2*)\\
\end{gathered}$$
Steady-state no volumetric source terms version
$$\begin{gathered}
{\left( {\rho {{\hat C}_p}} \right)_s}\left( {{\mathbf{v}}_s \cdot \nabla {T_s}} \right)- (1 - \varepsilon )\nabla \cdot {k_s}\nabla {T_s} - h({T_f} - {T_s}) = 0 \qquad (1) \\
+ {\left( {\rho {{\hat C}_p}} \right)_f}\left( {{\mathbf{v}} \cdot \nabla {T_f}} \right) - \varepsilon \nabla \cdot {k_f}\nabla {T_f} - h({T_s} - {T_f}) = 0 \qquad (2)\\
\end{gathered}$$
Normally, we assume that the solid is fixed $\left(v_s={0,0}\right)$, but in this case we will retain it to use it to create an interface mixing rule. Otherwise, the solid density and heat capacity drops out of the equation.
### Dimensional Analysis
When possible, I prefer to work with dimensionless equations. To make the system dimensionless, we divide the variables by characteristic variables and perform some factoring so the left hand and right hand side of the equations are dimensionless. Dimensionless temperatures and coordinates are defined as:
$$\begin{gathered}
{\Theta _s} = \frac{{{T_s} - {T_0}}}{{{T_c}}};{T_s} = {T_c}{\Theta _s} + T \\
{\Theta _f} = \frac{{{T_f} - {T_0}}}{{{T_c}}};{T_f} = {T_c}{\Theta _f} + T \\
xs = \frac{x}{D};x = D \cdot xs \\
ys = \frac{y}{D};y = D \cdot ys \\
T_c = \frac{q_{0}D}{k_f}
\end{gathered} $$
Making the appropriate substitutions, equations (1) and (2) become:
$$\begin{gathered}
\frac{{\left( {\rho {{\hat C}_p}} \right)_s}}{{\left( {\rho {{\hat C}_p}} \right)_f}}{{\mathbf{v}}^*}_s \cdot {\nabla ^*}{\Theta _s}- \left( {1 - \varepsilon } \right)\frac{{{k_s}}}{{{k_f}}}\frac{{{\tau _{conv}}}}{{{\tau _{cond}}}}{\nabla ^*}^2{\Theta _s} - \frac{{{\tau _{conv}}}}{{{\tau _{flux}}}}({\Theta _f} - {\Theta _s}) = 0 \qquad (3) \\
{{\mathbf{v}}^*} \cdot {\nabla ^*}{\Theta _f} - \varepsilon \frac{{{\tau _{conv}}}}{{{\tau _{cond}}}}{\nabla ^*}^2{\Theta _f} - \frac{{{\tau _{conv}}}}{{{\tau _{flux}}}}({\Theta _s} - {\Theta _f}) = 0 \qquad (4) \\
\end{gathered} $$
Where the various timescales are defined by:
$$\begin{gathered}
{\tau _{conv}} = \frac{D}{{{v_0}}} \\
{\tau _{cond}} = \frac{{{D^2}}}{\alpha };\alpha = \frac{{{k_f}}}{{{{\left( {\rho {{\hat C}_p}} \right)}_f}}} \\
{\tau _{flux}} = \frac{{{{\left( {\rho {{\hat C}_p}} \right)}_f}}}{h} \\
\end{gathered} $$
Consolidating the timescale ratios and noting that porosity and velocity are a function of the mesh region we are in.
$$\begin{gathered}
\kappa{{\mathbf{v}}^*}_s(\Omega ) \cdot {\nabla ^*}{\Theta _s}- \left( {1 - \varepsilon (\Omega )} \right)\sigma \gamma {\nabla ^*}^2{\Theta _s} - \eta ({\Theta _f} - {\Theta _s}) = 0 \qquad (5) \\
{{\mathbf{v}}^*}(\Omega ) \cdot {\nabla ^*}{\Theta _f} - \varepsilon (\Omega )\gamma {\nabla ^*}^2{\Theta _f} - \eta ({\Theta _s} - {\Theta _f}) = 0 \qquad (6) \\
\end{gathered} $$
Where
$$\begin{gathered}
\kappa = \frac{{\left( {\rho {{\hat C}_p}} \right)_s}}{{\left( {\rho {{\hat C}_p}} \right)_f}} \\
\sigma = \frac{{{k_s}}}{{{k_f}}} \\
\gamma = \frac{{{\tau _{conv}}}}{{{\tau _{cond}}}} \\
\eta = \frac{{{\tau _{conv}}}}{{{\tau _{flux}}}} \\
{{\mathbf{v}}^*} = \frac{{\mathbf{v}}}{{{v_0}}} \\
{{\mathbf{v}_s}^*} = \frac{{\mathbf{v}_s}}{{{v_0}}}
\end{gathered} $$
Now, we have a system of PDEs that depend on 5 dimensionless parameters.
### Neumann Boundary Condition (Flux) on Porous Boundary
When applying a flux on the porous boundary, there are different possible pathways through the fluid and solid. We must decide how to partition the heat flux. For parallel pathways, the effective conductivity is given by:
$${k_{eff}} = \left( {1 - \varepsilon } \right){k_s} + \varepsilon {k_f}$$
$$\frac{{{k_{eff}}}}{{{k_f}}} = \left( {1 - \varepsilon } \right)\sigma + \varepsilon $$
We will use this as a guide to partition the flux accordingly.
$$\begin{gathered}
{x_f} = \frac{\varepsilon }{{\left( {1 - \varepsilon } \right)\sigma + \varepsilon }} \\
{q_f} = {x_f}{q_{tot}} \\
{q_s} = (1 - {x_f}){q_{tot}} \\
\end{gathered} $$
## Modeling in _Mathematica_
With modeling, I am a firm believer that one should start simple, verify and validate the simulation results, and then build complexity. We will assume the convective diffusive heat equation has been previously validated.
Let's start with a model of the porous domain only.
### Mesh Creation
I followed the tutorial in the documentation for __[Element Mesh Creation](https://reference.wolfram.com/language/FEMDocumentation/tutorial/ElementMeshCreation.html)__ to generate a mesh tailored to the physical situation. You will need to load the FEM package like so
Needs["NDSolve`FEM`"]
Let's start with some definitions of a domain the is 1 unit high and 2 units long. We will also create some associations to collect surfaces and regions for clarity.
ht = 1;
len = 2;
top = ht;
bot = 0;
left = 0;
right = len;
bounds = <|inlet -> 1, hot -> 2|>;
regs = <|fluid -> 10, porous -> 20 , interface -> 15|>;
We can create a boundary mesh using the ToBoundaryMesh\[\] command and mark unique edges to be used as boundary conditions later. All other edges will be assigned the default boundary condtion of an insulated wall. The left edge colored $\color{Green}{Green}$ will describe the inlet condition and the bottom edge colored $\color{Red}{Red}$ will describe the heat flux condtion. All other boundaries will be the default Neumann zero flux condition.
bmesh = ToBoundaryMesh[
"Coordinates" -> {{left, bot}(*1*), {right, bot}(*2*), {right,
top}(*3*), {left, top}(*4*)},
"BoundaryElements" -> {LineElement[{{1, 2}(*bottom edge*)(*1*), {4,
1}(*left edge*)(*2*), {2, 3}(*3*), {3, 4}(*4*)}, {bounds[
hot], bounds[inlet], 3, 3}]}];
bmesh["Wireframe"["MeshElementMarkerStyle" -> Blue,
"MeshElementStyle" -> {Green, Red, Black}, ImageSize -> Large]]
#### Boundary Mesh Highlighting Special Boundaries
![enter image description here][2]
Now, we will create an element mesh that will have a $\color{Red}{Red}$ porous region with marker ID 20 (regs\[porous\]).
mesh = ToElementMesh[bmesh,
"RegionMarker" -> {{{(left + right)/2, (bot + top)/2},
regs[porous], 0.002}}];
mesh["Wireframe"["MeshElementStyle" -> {FaceForm[Red]},
ImageSize -> Large]]
#### Surface Mesh
![enter image description here][3]
# _Mathematica_ Implementation
By investing the time we did in mesh development, we can take advantage of some facilities provided by the Wolfram developers. First, we will convert our symbols to Latin for easier reading and define some initial parameters for testing.
$$\begin{gathered}
\varepsilon = eps = 0.4 \\
\sigma = s = 10 \\
\gamma = g = 2 \\
\eta = h = 0.2 \\
q_{tot} = 10 \\
v_{f} = \binom{6}{0} \\
v_{s} = \binom{0}{0}
\end{gathered} $$
eps = 0.4;
s = 10;
g = 2;
h = 1/5;
mu = 5000/h;
qtot = 10;
(* Derived Parameters *)
xf = eps/((1 - eps) s + eps);
hcapfrac = 5000;
qf = xf qtot;
qs = ( 1 - xf) qtot;
p = eps;
v = {6, 0};
vs = {0, 0};
fac = 1;
We can also set our Neumann/flux boundary conditions for each state variable on the hot wall.
nvsolid = NeumannValue[qs, ElementMarker == bounds[hot]];
nvfluid = NeumannValue[qf, ElementMarker == bounds[hot]];
Likewise, we can set the inlet boundary for each state variable to 0 with Dirichlet conditions
dcsolid =
DirichletCondition[ts[x, y] == 0, ElementMarker == bounds[inlet]];
dcfluid =
DirichletCondition[tf[x, y] == 0, ElementMarker == bounds[inlet]];
Recasting (5) and (6) into _Mathematica_ format (remember to put Neumann values on the right side of the equation where 0 was) we obtain:
fluideqn =
v.Inactive[Grad][tf[x, y], {x, y}] -
p g Inactive[Laplacian][tf[x, y], {x, y}] -
fac h (ts[x, y] - tf[x, y]) == nvfluid;
solideqn =
hcapfrac vs.Inactive[Grad][ts[x, y], {x, y}] - (1 -
p) g s Inactive[Laplacian][ts[x, y], {x, y}] -
fac h (tf[x, y] - ts[x, y]) == nvsolid;
Now, we can solve on the element mesh we created earlier and do a 3D plot of the solution (Solid Temperature is translucent and has dashed lines).
ifun = NDSolveValue[{fluideqn, solideqn, dcfluid, dcsolid}, {tf,
ts}, {x, y} \[Element] mesh];
opac = If[
ifun[[1]][0.5, 0.5] > ifun[[2]][0.5, 0.5], {0.25, 1}, {1, 0.25}];
pltf = Plot3D[ifun[[1]][x, y], {x, y} \[Element] mesh,
PlotRange -> All, PlotPoints -> {200, 200},
ColorFunction -> (Directive[Opacity[opac[[1]]],
ColorData["DarkBands"][#3]] &), MeshFunctions -> {#3 &},
Mesh -> 10, AxesLabel -> Automatic, MeshStyle -> {Black, Thick},
ImageSize -> Large];
plts = Plot3D[ifun[[2]][x, y], {x, y} \[Element] mesh,
PlotRange -> All, PlotPoints -> {200, 200},
ColorFunction -> (Directive[Opacity[opac[[2]]],
ColorData["DarkBands"][#3]] &), MeshFunctions -> {#3 &},
Mesh -> 10, AxesLabel -> Automatic,
MeshStyle -> {Black, Thick, Dashed}, ImageSize -> Large];
Grid[{{Show[{pltf, plts}, ViewProjection -> "Orthographic",
ViewPoint -> Front, ImageSize -> 450,
Background -> RGBColor[0.84`, 0.92`, 1.`], Boxed -> False],
Show[{pltf, plts}, ViewProjection -> "Orthographic",
ViewPoint -> Left , ImageSize -> 450,
Background -> RGBColor[0.84`, 0.92`, 1.`],
Boxed -> False] }, {Show[{pltf, plts},
ViewProjection -> "Orthographic", ViewPoint -> Top ,
ImageSize -> 450, Background -> RGBColor[0.84`, 0.92`, 1.`],
Boxed -> False],
Show[{pltf, plts}, ViewProjection -> "Perspective",
ViewPoint -> { Above, Right, Front} , ImageSize -> 450,
Background -> RGBColor[0.84`, 0.92`, 1.`], Boxed -> False]}},
Dividers -> Center]
#### Solution from Several Views
![enter image description here][4]
## Comparison to Another Code
Because we specified a low value of $\eta$, the heat transfer between the solid and fluid phase is poor leading to a large separation of temperatures. I think the plots look nice, but do they have any basis in reality? It is always good practice to validate models versus experiments or other codes when possible.
Fortunately, I have access to other codes and can say _Mathematica_ compares favorably as shown below ($T_{s} = Thick\ Lines$). Having such agreement gives me confidence to continue to the more complex model.
#### Results from Another Code for Porous Only Domain
![enter image description here][5]
We shoud note that ${\left. T_s\neq T_f \right|_{y = top}}$ at the top of the insulated wall.
# Porous Media + Fluid Domain
How energy is transferred through the interface between two types of media can be tricky business. Also, I don't believe that _Mathematica_ can have internal boundary conditions currently. However, _Mathematica_ provides other machinery to create an "interface". Recall, that I claimed that physically there must be a solid impermeable layer between the porous media and the fluid interface otherwise the porous media flow would short circuit into the fluid domain. So, let's include a thin layer in the model. This thin layer will allow us to equilibrate the solid and fluid temperatures by drastically increasing the fluid-solid heat transfer coefficient to the correct heat capacity weighted value and ramping the porosity to unity at the fluid interface so that we don't need to add another equation. A conceptual sketch is shown below.
#### Conceptual Model Zoomed into Interface Region
![enter image description here][6]
$$\begin{gathered}
\kappa{{\mathbf{v}}^*}_s(\Omega ) \cdot {\nabla ^*}{\Theta _s}- \left( {1 - \varepsilon (\Omega )} \right)\sigma \gamma {\nabla ^*}^2{\Theta _s} - -\mu(\Omega )\eta ({\Theta _f} - {\Theta _s}) = 0 \qquad (7) \\
{{\mathbf{v}}^*}(\Omega ) \cdot {\nabla ^*}{\Theta _f} - \varepsilon (\Omega )\gamma {\nabla ^*}^2{\Theta _f} -\mu(\Omega ) \eta ({\Theta _s} - {\Theta _f}) = 0 \qquad (8) \\
\end{gathered} $$
Where
$${{\mathbf{v}}^*}_s(\Omega )=
\begin{Bmatrix}
\binom{0}{0.001} & \Omega=Interface\\
\binom{0}{0} & True
\end{Bmatrix}$$
$${{\mathbf{v}}^*}(\Omega )=
\begin{Bmatrix}
\binom{6}{0} & \Omega=Porous\\
\binom{20 (1 - (\frac{y}{r})^2}{0} & \Omega=Fluid\\
\binom{0}{0.001} & \Omega=Interface\\
\end{Bmatrix}$$
$$\varepsilon(\Omega )=
\begin{Bmatrix}
\varepsilon_0 & \Omega=Porous\\
\frac{\left(1-\varepsilon_0\right) (r+thick+y)}{thick}+\varepsilon _0 & \Omega=Interface\\
1 & True
\end{Bmatrix}$$
$$\mu(\Omega )=
\begin{Bmatrix}
\frac{500,000}{\eta} & \Omega=Interface\\
1 & True
\end{Bmatrix}$$
### Mesh Creation for Porous, Fluid, and Coupling Interface
We will create a domain that is 1 unit in height by 2 units in length. The fluid ($\color{Green}{Green}$) will occupy the top 1/4 of the domain. The interface layer ($\color{Cyan}{Cyan}$) will be 1/40 the height. We will mark special boundaries and regions with unique IDs for boundary condition assignment and to change properties and mesh refine regions.
ht = 1;
len = 2;
r = ht/8;
top = r;
bot = r - ht;
left = 0;
right = len;
thick = ht/40;
interfacef = -r;
interfaces = -r - thick;
buffer = 1.5 thick;
(* Use associations for clearer assignment later *)
bounds = <|inlet -> 1, hot -> 2|>;
regs = <|porous -> 10, fluid -> 20, interface -> 15|>;
(* Meshing Definitions *)
(* Coordinates *)
crds = {{left, bot}(*1*), {right, bot}(*2*), {right, top}(*3*), {left,
top}(*4*), {left, interfacef}(*5*), {right,
interfacef}(*6*), {left, interfaces}(*7*), {right, interfaces}(*8*)};
(* Edges *)
lelms = {{1, 7}, {7, 5}, {5, 4}(*left edge*)(*1*), {1,
2}(*bottom edge*)(*2*), {2, 8}, {8, 6}, {6, 3}, {3, 4}, {5,
6}, {7, 8}(*3*)};
boundaryMarker = {bounds[inlet], bounds[inlet], bounds[inlet],
bounds[hot], 3, 3, 3, 3, 3, 3}; (* 3 will be a default boundary *)
bcEle = {LineElement[lelms, boundaryMarker]};
bmesh = ToBoundaryMesh["Coordinates" -> crds,
"BoundaryElements" -> bcEle];
bmesh["Wireframe"["MeshElementMarkerStyle" -> Blue,
"MeshElementStyle" -> {Green, Red, Black}, ImageSize -> Large]]
(* 2D Regions *)
(* Identify Center Points of Different Material Regions *)
fluidCenter = {(left + right)/2, 0};
fluidReg = {fluidCenter, regs[fluid], 0.0005};
interfaceCenter = {(left + right)/2, (interfaces + interfacef)/2};
interfaceReg = {interfaceCenter, regs[interface], 0.00005};
porousCenter = {(left + right)/2, (bot + interfaces)/2};
porousReg = {porousCenter, regs[porous], 0.002};
meshRegs = {fluidReg, interfaceReg, porousReg};
mesh = ToElementMesh[bmesh, "RegionMarker" -> meshRegs,
MeshRefinementFunction ->
Function[{vertices, area}, Block[{x, y}, {x, y} = Mean[vertices];
If[y > (interfaces + interfacef)/2 - buffer &&
y < (interfaces + interfacef)/2 + buffer, area > 0.00005,
area > 0.01]]]];
mesh["Wireframe"[
"MeshElementStyle" -> {FaceForm[Green], FaceForm[Yellow],
FaceForm[Red]}, ImageSize -> Large]]
Show[mesh[
"Wireframe"[
"MeshElementStyle" -> {FaceForm[Green], FaceForm[Yellow],
FaceForm[Red]}, ImageSize -> Large]],
PlotRange -> {{0,
0.5}, {(interfaces + interfacef)/2 -
4 buffer, (interfaces + interfacef)/2 + 4 buffer}}]
#### Boundary Mesh
![enter image description here][7]
#### Full Mesh
![enter image description here][8]
#### Mesh Zoomed into Interface Region
![enter image description here][9]
### NDSolve Setup and Solution
eps = 0.4;
s = 10;
g = 2;
h = 1/5;
mu = 500000/h;
qtot = 10;
xf = eps/((1 - eps) s + eps);
kappa = 500;
qf = xf qtot;
qs = ( 1 - xf) qtot;
(* Region Dependent Properties with Piecewise Functions *)
p = Evaluate[Piecewise[{{eps, ElementMarker == regs[porous]},
{(1 - eps)/thick (y + (r + thick)) + eps,
ElementMarker == regs[interface]},
{1, True}}]];
v = Evaluate[Piecewise[{{{6, 0}, ElementMarker == regs[porous]},
{{20 (1 - (y/r)^2), 0}, ElementMarker == regs[fluid]},
{{0, 0.0001}, ElementMarker == regs[interface]}}]];
vs = Evaluate[
Piecewise[{{{0, 0.0001}, ElementMarker == regs[interface]},
{{0, 0}, True}}]];
(* fac increases heat transfer coefficient in interface layer *)
fac = Evaluate[If[ElementMarker == regs[interface], mu, 1]];
(* Neumann Conditions for the Hot Bottom Wall *)
nvsolid = NeumannValue[qs, ElementMarker == bounds[hot]];
nvfluid = NeumannValue[qf, ElementMarker == bounds[hot]];
(* Dirichlet Conditions for the Left Wall *)
dcsolid =
DirichletCondition[ts[x, y] == 0, ElementMarker == bounds[inlet]];
dcfluid =
DirichletCondition[tf[x, y] == 0, ElementMarker == bounds[inlet]];
(* Balance Equations for Fluid and Solid Temperature *)
fluideqn =
v.Inactive[Grad][tf[x, y], {x, y}] -
p g Inactive[Laplacian][tf[x, y], {x, y}] -
fac h (ts[x, y] - tf[x, y]) == nvfluid;
solideqn =
kappa vs.Inactive[Grad][ts[x, y], {x, y}] - (1 - p) g s Inactive[
Laplacian][ts[x, y], {x, y}] - fac h (tf[x, y] - ts[x, y]) ==
nvsolid;
ifun = NDSolveValue[{fluideqn, solideqn, dcfluid, dcsolid}, {tf,
ts}, {x, y} \[Element] mesh];
(* Set opacity for higher temperature to be translucent so we can see \
the layer below *)
opac = If[
ifun[[1]][0.5, (-r + bot)/2] > ifun[[2]][0.5, (-r + bot)/2], {0.25,
1}, {1, 0.25}];
(* Plots *)
pltf = Plot3D[ifun[[1]][x, y], {x, y} \[Element] mesh,
PlotRange -> All, PlotPoints -> {200, 200},
ColorFunction -> (Directive[Opacity[opac[[1]]],
ColorData["DarkBands"][#3]] &), MeshFunctions -> {#3 &},
Mesh -> 18, AxesLabel -> Automatic, MeshStyle -> {Black, Thick},
ImageSize -> Large];
plts = Plot3D[ifun[[2]][x, y], {x, y} \[Element] mesh,
PlotRange -> All, PlotPoints -> {200, 200},
ColorFunction -> (Directive[Opacity[opac[[2]]],
ColorData["DarkBands"][#3]] &), MeshFunctions -> {#3 &},
Mesh -> 18, AxesLabel -> Automatic,
MeshStyle -> {Black, Thick, Dashed}, ImageSize -> Large];
Grid[{{Show[{pltf, plts}, ViewProjection -> "Orthographic",
ViewPoint -> Front, ImageSize -> 400,
Background -> RGBColor[0.84`, 0.92`, 1.`], Boxed -> False],
Show[{pltf, plts}, ViewProjection -> "Orthographic",
ViewPoint -> Left , ImageSize -> 400,
Background -> RGBColor[0.84`, 0.92`, 1.`],
Boxed -> False] }, {Show[{pltf, plts},
ViewProjection -> "Orthographic", ViewPoint -> Top ,
ImageSize -> 400, Background -> RGBColor[0.84`, 0.92`, 1.`],
Boxed -> False],
Show[{pltf, plts}, ViewProjection -> "Perspective",
ViewPoint -> { Above, Right, Front} , ImageSize -> 400,
Background -> RGBColor[0.84`, 0.92`, 1.`], Boxed -> False]}},
Dividers -> Center]
#### Full Model Porous + Interface + Fluid
![enter image description here][10]
## Comparison to Another Code
Again, as shown below, the _Mathematica_ solution compares qualitatively favorably to another code. The contour plots have similar shapes but _Mathematica_ has about 40% higher temperatures, which could indicate that I dropped or added a $\varepsilon$ somewhere in the setup of the other code, but a detailed debugging is probably in order. The implementation in the other code is slightly different due to different capabilities which I will elaborate. In the other code, the interface layer is a just a high thermal conductivity solid. I am able to set the porous-solid interface temperature to be equal to the porous temperature (solid-fluid bulk mix temperature).
#### Another Code Full Model Porous + Interface + Fluid
![enter image description here][11]
Finally, because I know that I can set an interface boundary condition in the other code, I can remove the interface region and set the fluid boundary condition to the porous mean temperature and judge the effect. As shown in the image below, it probably would not change conclusions much.
#### Another Code Full Model Porous + Fluid No Interface
![enter image description here][12]
# Future Enhancements
After looking at the documentation and reading the literature, I should recast the equations to conform to the coefficient form of the PDE so that I could additionally scale the $x$ coordinate by the Péclet number. This would allow for easier comparison to literature and keep a more reasonable aspect ratio for the domain.
$$\nabla \cdot\left( { - \mathbf{c}\nabla u - \alpha u + \gamma } \right) + \beta \cdot\nabla u + au - f = 0$$
# Summary
- An exact analytical solution to an LTNE problem coupled to a fluid domain seems improbable.
- _Mathematica_
- Currently, I don't think interface BCs can be defined in the FEM framework.
- Does provide meshing features that allow the construction of an interface with finite thickness that might be used a workaround to this limitation.
- Results
- Porous Media Only Domain.
- Quite easy to setup a mesh and identify unique boundaries and regions.
- Straightforward to translate our system of PDEs created by hand into equations suitable for _Mathematica_.
- Relatively easy to post process results using examples from documentation.
- Solution compares favorably to another PDE code.
- Combined Porous Media--Interface--Fluid Domains
- Quite easy to setup a mesh with appropriate refinement and identify unique boundaries and regions.
- Straightforward to translate our system of PDEs created by hand into equations suitable for _Mathematica_.
- Relatively easy to post process results using examples from documentation.
- Solution compares qualitatively favorably to another PDE code.
- Additional study/debugging will be required to fully harmonize the results.
- Could be a simple matter of missing a parameter or something deeper.
- Conclusion
- A prototype of a simple LTNE model was developed in _Mathematica_.
- In modeling, there are many ways to skin a cat and this is just one possible way that produces a solution in a resonable period. More elegant solutions are welcome.
- Additional model development and validation is necessary.
- _Mathematica_ provides enough features and tools that you can get close to what you need with some imagination.
# Attachment
The attached notebook has the complete code to produce a solution to study and modify as one sees fit. I did recast the equations into coefficient form because that appears to be better practice.
Storing the solution will increase notebook size >100MB.
[1]: http://community.wolfram.com//c/portal/getImageAttachment?filename=SystemDescription.png&userId=1402928
[2]: http://community.wolfram.com//c/portal/getImageAttachment?filename=porousonlyboundaries.png&userId=1402928
[3]: http://community.wolfram.com//c/portal/getImageAttachment?filename=porousonlymesh.png&userId=1402928
[4]: http://community.wolfram.com//c/portal/getImageAttachment?filename=porousonlyPlot3D.png&userId=1402928
[5]: http://community.wolfram.com//c/portal/getImageAttachment?filename=LTNEPorousOnly2.png&userId=1402928
[6]: http://community.wolfram.com//c/portal/getImageAttachment?filename=InterfaceLayer.png&userId=1402928
[7]: http://community.wolfram.com//c/portal/getImageAttachment?filename=fullmodelboundaries.png&userId=1402928
[8]: http://community.wolfram.com//c/portal/getImageAttachment?filename=fullmodelmesh.png&userId=1402928
[9]: http://community.wolfram.com//c/portal/getImageAttachment?filename=fullmodelzoom.png&userId=1402928
[10]: http://community.wolfram.com//c/portal/getImageAttachment?filename=PorousMediaSimulation.png&userId=1402928
[11]: http://community.wolfram.com//c/portal/getImageAttachment?filename=LTNEPConductorFluid4.png&userId=1402928
[12]: http://community.wolfram.com//c/portal/getImageAttachment?filename=LTNEPorousFluidNoInterface.png&userId=1402928Tim Laska2018-09-20T17:43:10ZWhich Riemann Zeta function representation is used in Zeta[ ] ?
http://community.wolfram.com/groups/-/m/t/1467997
Which representation of the Riemann Zeta function does Mathematica use for zeta[ ] ?Rick Jarosh2018-09-20T02:26:11ZThoughts on a Python interface, and why ExternalEvaluate is just not enough
http://community.wolfram.com/groups/-/m/t/1185247
`ExternalEvaluate`, introduced in M11.2, is a nice initiative. It enables limited communication with multiple languages, including Python, and appears to be designed to be relatively easily extensible (see ``ExternalEvaluate`AddHeuristic`` if you want to investigate, though I wouldn't invest in this until it becomes documented).
**My great fear, however, is that with `ExternalEvaluate` Wolfram will consider the question of a Python interface settled.**
This would be a big mistake. A *general* framework, like `ExternalEvaluate`, that aims to work with *any* language and relies on passing code (contained in a string) to an evaluator and getting JSON back, will never be fast enough or flexible enough for *practical scientific computing*.
Consider a task as simple as computing the inverse of a $100\times100$ Mathematica matrix using Python (using [`numpy.linalg.inv`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.inv.html)).
I challenge people to implement this with `ExternalEvaluate`. It's not possible to do it *in a practically useful way*. The matrix has to be sent *as code*, and piecing together code from strings just can't replace structured communication. The result will need to be received as something encodable to JSON. This has terrible performance due to multiple conversions, and even risks losing numerical precision.
Just sending and receiving a tiny list of 10000 integers takes half a second (!)
In[6]:= ExternalEvaluate[py, "range(10000)"]; // AbsoluteTiming
Out[6]= {0.52292, Null}
Since I am primarily interested in scientific and numerical computing (as I believe most M users are), I simply won't use `ExternalEvaluate` much, as it's not suitable for this purpose. What if we need to do a [mesh transformation](https://mathematica.stackexchange.com/q/155484/12) that Mathematica can't currently handle, but there's a Python package for it? It's exactly the kind of problem I am looking to apply Python for. I have in fact done mesh transformations using MATLAB toolboxes directly from within Mathematica, using [MATLink][1], while doing the rest of the processing in Mathematica. But I couldn't do this with ExternalEvaluate/Python in a reasonable way.
In 2017, any scientific computing system *needs* to have a Python interface to be taken seriously. [MATLAB has one][2], and it *is* practically usable for numerical/scientific problems.
----
## A Python interface
I envision a Python interface which works like this:
- The MathLink/WSTP API is exposed to Python, and serves as the basis of the system. MathLink is good at transferring large numerical arrays efficiently.
- Fundamental data types (lists, dictionaries, bignums, etc.) as well as datatypes critical for numerical computing (numpy arrays) can be transferred *efficiently* and *bidirectionally*. Numpy arrays in particular must translate to/from packed arrays in Mathematica with the lowest possible overhead.
- Python functions can be set up to be called from within Mathematica, with automatic argument translation and return type translation. E.g.,
PyFun["myfun"][ (* myfun is a function defined in Python *)
{1,2,3} (* a list *),
PyNum[{1,2,3}] (* cast to numpy array, since the interpretation of {1,2,3} is ambiguous *),
PySet[{1,2,3}] (* cast to a set *)
]
- The system should be user-extensible to add translations for new datatypes, e.g. a Python class that is needed frequently for some application.
- The primary mode of operation should be that Python is run as a slave (subprocess) of Mathematica. But there should be a second mode of operation where both Mathematica and Python are being used interactively, and they are able to send/receive structured data to/from each other on demand.
- As a bonus: Python can also call back to Mathematica, so e.g. we can use a numerical optimizer available in Python to find the minimum of a function defined in Mathematica
- An interface whose primary purpose is to call Mathematica from Python is a different topic, but can be built on the same data translation framework described above.
The development of such an interface should be driven by real use cases. Ideally, Wolfram should talk to users who use Mathematica for more than fun and games, and do scientific computing as part of their daily work, with multiple tools (not just M). Start with a number of realistic problems, and make sure the interface can help in solving them. As a non-trivial test case for the datatype-extension framework, make sure people can set up auto-translation for [SymPy objects][3], or a [Pandas dataframe][4], or a [networkx graph][5]. Run `FindMinimum` on a Python function and make sure it performs well. (In a practical scenario this could be a function implementing a physics simulation rather than a simple formula.) As a performance stress test, run `Plot3D` (which triggers a very high number of evaluations) on a Python function. Performance and usability problems will be exposed by such testing early, and then the interface can be *designed* in such a way as to make these problems at least solvable (if not immediately solved in the first version). I do not believe that they are solvable with the `ExternalEvaluate` design.
Of course, this is not the only possible design for an interface. J/Link works differently: it has handles to Java-side objects. But it also has a different goal. Based on my experience with MATLink and RLink, I believe that *for practical scientific/numerical computing*, the right approach is what I outlined above, and that the performance of data structre translation is critical.
----
## ExternalEvaluate
Don't get me wrong, I do think that the `ExternalEvaluate` framework is a very useful initiative, and it has its place. I am saying this because I looked at its source code and it appears to be easily extensible. R has zeromq and JSON capabilities, and it looks like one could set it up to work with `ExternalEvaluate` in a day or so. So does Perl, anyone want to give it a try? `ExternalEvaluate` is great because it is simple to use and works (or can be made to work) with just about any interpreted language that speaks JSON and zeromq. But it is also, in essence, a quick and dirty hack (that's extensible in a quick and dirty way), and won't be able to scale to the types of problems I mentioned above.
----
## MathLink/WSTP
Let me finally say a few words about why MathLink/WSTP are critical for Mathematica, and what should be improved about them.
I believe that any serious interface should be built on top of MathLink. Since Mathematica already has a good interface capable of inter-process communication, that is designed to work well with Mathematica, and designed to handle numerical and symbolic data efficiently, use it!!
Two things are missing:
- Better documentation and example programs, so more people will learn MathLink
- If the MathLink library (not Mathematica!) were open source, people would be able to use it to link to libraries [which are licensed under the GPL][6]. Even a separate open source implementation that only supports shared memory passing would be sufficient—no need to publish the currently used code in full. Many scientific libraries are licensed under the GPL, often without their authors even realizing that they are practically preventing them from being used from closed source systems like Mathematica (due to the need to link to the MathLink libraries). To be precise, GPL licensed code *can* be linked with Mathematica, but the result cannot be shared with anyone. I have personally requested the author of a certain library to grant an exception for linking to Mathematica, and they did not grant it. Even worse, I am not sure they understood the issue. The authors of other libraries *cannot* grant such a permission because they themselves are using yet other GPL's libraries.
[MathLink already has a more permissive license than Mathematica.][7] Why not go all the way and publish an open source implementation?
I am hoping that Wolfram will fix these two problems, and encourage people to create MathLink-based interfaces to other systems. (However, I also hope that Wolfram will create a high-quality Python link themselves instead of relying on the community.)
I have talked about the potential of Mathematica as a glue-language at some Wolfram events in France, and I believe that the capability to interface external libraries/systems easily is critical for Mathematica's future, and so is a healthy third-party package ecosystem.
[1]: http://matlink.org/
[2]: https://www.mathworks.com/help/matlab/matlab-engine-for-python.html
[3]: http://www.sympy.org/
[4]: http://pandas.pydata.org/
[5]: https://networkx.github.io/
[6]: https://en.wikipedia.org/wiki/Copyleft
[7]: https://www.wolfram.com/legal/agreements/mathlink.htmlSzabolcs Horvát2017-09-15T12:33:04ZCan one have custom "operator records" in WSM?
http://community.wolfram.com/groups/-/m/t/1455024
Inspired by [this answer][1] to my question on stack overflow I was curious about overloading the class **operator** to have an **operator record**. From the [Modelica specs for Version 3.2.2][2] (Chapter 14) I take it that *operator records* are available in the language since 2013.
**Can one have a class *operator* in WSM?**
Edit:
I edited the title to make my question more precise: I am interested in defining my own operator records.
Update: I have cross-posted this question now on Mathematica Stack Exchange [(181939)][3].
[1]: https://stackoverflow.com/a/52313716/5363743
[2]: https://www.modelica.org/documents/ModelicaSpec32Revision2.pdf
[3]: https://mathematica.stackexchange.com/q/181939/764Guido Wolf Reichert2018-09-13T14:27:01ZMetaprogramming: the Future of the Wolfram Language
http://community.wolfram.com/groups/-/m/t/1435093
With all the marvelous new functionality that we have come to expect with each release, it is sometimes challenging to maintain a grasp on what the Wolfram language encompasses currently, let alone imagine what it might look like in another ten years. Indeed, the pace of development appears to be accelerating, rather than slowing down.
However, I predict that the "problem" is soon about to get much, much worse. What I foresee is a step change in the pace of development of the Wolfram Language that will produce in days and weeks, or perhaps even hours and minutes, functionality might currently take months or years to develop.
So obvious and clear cut is this development that I have hesitated to write about it, concerned that I am simply stating something that is blindingly obvious to everyone. But I have yet to see it even hinted at by others, including Wolfram. I find this surprising, because it will revolutionize the way in which not only the Wolfram language is developed in future, but in all likelihood programming and language development in general.
The key to this paradigm shift lies in the following unremarkable-looking WL function WolframLanguageData[], which gives a list of all Wolfram Language symbols and their properties. So, for example, we have:
WolframLanguageData["SampleEntities"]
![enter image description here][1]
This means we can treat WL language constructs as objects, query their properties and apply functions to them, such as, for example:
WolframLanguageData["Cos", "RelationshipCommunityGraph"]
![enter image description here][2]
In other words, the WL gives us the ability to traverse the entirety of the WL itself, combining WL objects into expressions, or programs. This process is one definition of the term “Metaprogramming”.
What I am suggesting is that in future much of the heavy lifting will be carried out, not by developers, but by WL programs designed to produce code by metaprogramming. If successful, such an approach could streamline and accelerate the development process, speeding it up many times and, eventually, opening up areas of development that are currently beyond our imagination (and, possibly, our comprehension).
So how does one build a metaprogramming system? This is where I should hand off to a computer scientist (and will happily do so as soon as one steps forward to take up the discussion). But here is a simple outline of one approach.
The principal tool one might use for such a task is genetic programming:
WikipediaData["Genetic Programming"]
> In artificial intelligence, genetic programming (GP) is a technique whereby computer programs are encoded as a set of genes that are then modified (evolved) using an evolutionary algorithm (often a genetic algorithm, "GA") – it is an application of (for example) genetic algorithms where the space of solutions consists of computer programs. The results are computer programs that are able to perform well in a predefined task. The methods used to encode a computer program in an artificial chromosome and to evaluate its fitness with respect to the predefined task are central in the GP technique and still the subject of active research.
One can take issue with this explanation on several fronts, in particular the suggestion that GP is used primarily as a means of generating a computer program for performing a predefined task. That may certainly be the case, but need not be.
Leaving that aside, the idea in simple terms is that we write a program that traverses the WL structure in some way, splicing together language objects to create a WL program that “does something”. That “something” may be a predefined task and indeed this would be a great place to start: to write a GP metaprogramming system that creates WL programs that replicate the functionality of existing WL functions. Most of the generated programs would likely be uninteresting, slower versions of existing functions; but it is conceivable, I suppose, that some of the results might be of academic interest, or indicate a potentially faster computation method, perhaps. However, the point of the exercise is to get started on the metaprogramming project, with a simple(ish) task with very clear, pre-defined goals and producing results that are easily tested. In this case the “objective function” is a comparison of results produced by the inbuilt WL functions vs the GP-generated functions, across some selected domain for the inputs.
I glossed over the question of exactly how one “traverses the WL structure” for good reason: I feel sure that there must have been tremendous advances in the theory of how to do this in the last 50 years. But just to get the ball rolling, one could, for instance, operate a dual search, with a local search evaluating all of the functions closely connected to the (randomly chosen) starting function (WL object), while a second “long distance” search jumps randomly to a group of functions some specified number of steps away from the starting function.
[At this point I envisage the computer scientists rolling their eyes and muttering “doesn’t this idiot know about the {fill in the bank} theorem about efficient domain search algorithms?”].
Anyway, to continue. The initial exercise is about the mechanics of the process rather that the outcome. The second stage is much more challenging, as the goal is to develop new functionality, rather than simply to replicate what already exists. It would entail defining a much more complex objective function, as well as perhaps some constraints on program size, the number and types of WL objects used, etc.
An interesting exercise, for example, would be to try to develop a metaprogramming system capable of winning the Wolfram One-Liner contest. Here, one might characterize the objective function as “something interesting and surprising”, and we would impose a tight constraint on the length of programs generated by the metaprogramming system to a single line of code.
What is “interesting and surprising”? To be defined – that’s a central part of the challenge. But, in principle, I suppose one might try to train a neural network to classify whether or not a result is “interesting” based on the results of prior one-liner competitions.
From there, it’s on to the hard stuff: designing metaprogramming systems to produce WL programs of arbitrary length and complexity to do “interesting stuff” in a specific domain. That “interesting stuff” could be, for instance, a more efficient approximation for a certain type of computation, a new algorithm for detecting certain patterns, or coming up with some completely novel formula or computational concept.
Obviously one faces huge challenges in this undertaking; but the potential rewards are also enormous in terms of accelerating the pace of language development and discovery. It is a fascinating area for R&D, one that the WL is ideally situated to exploit. Indeed, I would be mightily surprised to learn that there is not already a team engaged on just such research at Wolfram. If so, perhaps one of them could comment here?
[1]: http://community.wolfram.com//c/portal/getImageAttachment?filename=10942Fig1.png&userId=773999
[2]: http://community.wolfram.com//c/portal/getImageAttachment?filename=O_12.png&userId=773999Jonathan Kinlay2018-09-02T13:38:13ZDecouple the following equations?
http://community.wolfram.com/groups/-/m/t/1459261
I have two couple equations and I want to decouple it. How can I do it? any help will be appreciated.
eq1 = k*X''[Y] + Bi*(Z[Y] - X[Y]) == Subscript[U, p]/Subscript[U, m];
eq2 = Z''[Y] - Bi*(Z[Y] - X[Y]) == 0;Mirza Farrukh Baig2018-09-15T13:55:30ZExtract the real solutions among many which are complex?
http://community.wolfram.com/groups/-/m/t/1450000
Here is the equation and the solution: I want the two real solutions.
In[2]:= soln1[a_, b_, c_, \[Xi]_, \[Lambda]_, \[Phi]_, rn_] :=
NSolve[\[Rho]^2 +
2 \[Rho] (a Sin[\[Theta]] Cos[\[Phi]] +
b Sin[\[Theta]] Sin[\[Phi]] + c Cos[\[Theta]]) + a^2 + b^2 +
c^2 - rn^2 ==
0 && \[Rho] - \[Lambda] rn Cos[\[Theta] + \[Xi] Cos[\[Phi]]]^2 ==
0, {\[Rho], \[Theta]}];
In[8]:= ss = soln1[.5, .5, .5, .2, .5, .1, 1]
During evaluation of In[8]:= NSolve::ifun: Inverse functions are being used by NSolve, so some solutions may not be found; use Reduce for complete solution information.
Out[8]= {{\[Rho] -> -2.51366 - 3.61652 I, \[Theta] ->
1.8313 - 1.79756 I}, {\[Rho] -> -2.51366 + 3.61652 I, \[Theta] ->
1.8313 + 1.79756 I}, {\[Rho] -> -0.146207 -
0.0312241 I, \[Theta] -> -1.82025 +
0.520643 I}, {\[Rho] -> -0.146207 +
0.0312241 I, \[Theta] -> -1.82025 - 0.520643 I}, {\[Rho] ->
0.152974, \[Theta] -> 0.785684}, {\[Rho] ->
0.311045 - 0.210578 I, \[Theta] ->
2.25043 - 0.388346 I}, {\[Rho] ->
0.311045 + 0.210578 I, \[Theta] ->
2.25043 + 0.388346 I}, {\[Rho] -> 0.417435, \[Theta] -> -0.617471}}Hong-Yee Chiu2018-09-11T16:44:53ZSolve the following non-linear integer max-problem with NMaximize?
http://community.wolfram.com/groups/-/m/t/1463232
Hi!
I have tried to solve the following non-linear integer max-problem, without setting some lower (positive) bounds on Q1, Q2 & Q3, just being positive. Mathematica fails to find the optimal solotion, unless I increase the lower bounds for Q2 & Q3 to about 40. I have solved the same problem in LINGO and found the global solution {P -> 60., Q1 -> 0, Q2 -> 52., Q3 -> 60., y1 -> 0, y2 -> 1, y3 -> 1, \z1 -> 0, z2 -> 0, z3 -> 1} in less than a second. Any suggestions?
NMaximize[{(Q1 + Q2 + Q3) (P - 20) + 3*0.5 (80 - P) Q1*z1 +
2*0.5 (100 - 0.8 P) Q2*z2 + 2*0.5 (90 - 0.5 P) Q3*z3,
z1 + z2 + z3 == 1 && Q1 == (80 - P) y1 && Q2 == (100 - 0.8 P) y2 &&
Q3 == (90 - 0.5 P) y3 && 2 <= y1 + y2 + y3 <= 3 &&
3 <= y1 + y2 + y3 + z1 + z2 + z3 <= 4 &&
2 <= y1 + y2 + y3 + z1 <= 3 && 2 <= y1 + y2 + y3 + z2 <= 3 &&
2 <= y1 + y2 + y3 + z3 <= 3 && Q1 >= 0 && Q2 >= 0 && Q3 >= 0 &&
1 >= y1 >= 0 && 1 >= y2 >= 0 && 1 >= y3 >= 0 && z1 >= 0 &&
z2 >= 0 && z3 >= 0 && 80 > P > 20 && y1 \[Element] Integers &&
y2 \[Element] Integers && y3 \[Element] Integers &&
z1 \[Element] Integers && z2 \[Element] Integers &&
z3 \[Element] Integers}, {P, Q1, Q2, Q3, y1, y2, y3, z1, z2, z3}]Christos Papahristodoulou2018-09-17T19:22:37ZFind the root of a numerical function that contains a root-finding itself?
http://community.wolfram.com/groups/-/m/t/1463531
I have encountered the following problem:
- define a numerical function h with input x and output y
- use FindRoot to find the root x* of h
- one characteristic of h: at some part within h FindRoot is called
I have done that with other software, but its not working in Mathematica, or I am doing something wrong. I have attached a MWE.
Simply evaluating the function h works, but FindRoot does return many error messages. I really do not understand why.
Any ideas?
Best,
Benjamin
h[inp_] :=
Module[{y, x, a, equ, sol, z},
a = inp;
equ = {x + .5 y - a[[1]], x + y - a[[2]]};
sol = FindRoot[equ, {{x, 1.0}, {y, 1.0}}];
z = {x, y} /. sol;
z - {1, 1}
]
XY = {X, Y};
XYval = {1.5, 2.0};
XYStart = {{XY[[1]], XYval[[1]]}, {XY[[2]], XYval[[2]]}};
h[XYval]
FindRoot[h[XYInput], XYStart]Benjamin L2018-09-17T20:28:02ZSolver for unsteady flow with the use of Mathematica FEM
http://community.wolfram.com/groups/-/m/t/1433064
![fig7][331]
I started the discussion [here][1] but I also want to repeat on this forum.
There are many commercial and open code for solving the problems of unsteady flows.
We are interested in the possibility of solving these problems using Mathematica FEM. Previously proposed solvers for stationary incompressible isothermal flows:
Solving 2D Incompressible Flows using Finite Elements:
http://community.wolfram.com/groups/-/m/t/610335
FEM Solver for Navier-Stokes equations in 2D:
http://community.wolfram.com/groups/-/m/t/611304
Nonlinear FEM Solver for Navier-Stokes equations in 2D:
https://mathematica.stackexchange.com/questions/94914/nonlinear-fem-solver-for-navier-stokes-equations-in-2d/96579#96579
We give several examples of the successful application of the finite element method for solving unsteady problem including nonisothermal and compressible flows. We will begin with two standard tests that were proposed to solve this class of problems by
M. Schäfer and S. Turek, Benchmark computations of laminar ﬂow around a cylinder (With support by F. Durst, E. Krause and R. Rannacher). In E. Hirschel, editor, Flow Simulation with High-Performance Computers II. DFG priority research program results 1993-1995, number 52 in Notes Numer. Fluid Mech., pp.547–566. Vieweg, Weisbaden, 1996. https://www.uio.no/studier/emner/matnat/math/MEK4300/v14/undervisningsmateriale/schaeferturek1996.pdf
![fig8][332]
Let us consider the flow in a flat channel around a cylinder at Reynolds number = 100, when self-oscillations occur leading to the detachment of vortices in the aft part of cylinder. In this problem it is necessary to calculate drag coeﬃcient, lift coeﬃcient and pressure diﬀerence in the frontal and aft part of the cylinder as functions of time, maximum drag coeﬃcient, maximum lift coeﬃcient , Strouhal number and pressure diﬀerence $\Delta P(t)$ at $t = t0 +1/2f$. The frequency f is determined by the period of oscillations of lift coeﬃcient f=f(c_L). The data for this test, the code and the results are shown below.
H = .41; L = 2.2; {x0, y0, r0} = {1/5, 1/5, 1/20};
Ω = RegionDifference[Rectangle[{0, 0}, {L, H}], Disk[{x0, y0}, r0]];
RegionPlot[Ω, AspectRatio -> Automatic]
K = 2000; Um = 1.5; ν = 10^-3; t0 = .004;
U0[y_, t_] := 4*Um*y/H*(1 - y/H)
UX[0][x_, y_] := 0;
VY[0][x_, y_] := 0;
P0[0][x_, y_] := 0;
Do[
{UX[i], VY[i], P0[i]} =
NDSolveValue[{{Inactive[
Div][({{-μ, 0}, {0, -μ}}.Inactive[Grad][
u[x, y], {x, y}]), {x, y}] +
\!\(\*SuperscriptBox[\(p\),
TagBox[
RowBox[{"(",
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y] + (u[x, y] - UX[i - 1][x, y])/t0 +
UX[i - 1][x, y]*D[u[x, y], x] +
VY[i - 1][x, y]*D[u[x, y], y],
Inactive[
Div][({{-μ, 0}, {0, -μ}}.Inactive[Grad][
v[x, y], {x, y}]), {x, y}] +
\!\(\*SuperscriptBox[\(p\),
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y] + (v[x, y] - VY[i - 1][x, y])/t0 +
UX[i - 1][x, y]*D[v[x, y], x] +
VY[i - 1][x, y]*D[v[x, y], y],
\!\(\*SuperscriptBox[\(u\),
TagBox[
RowBox[{"(",
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y] +
\!\(\*SuperscriptBox[\(v\),
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y]} == {0, 0, 0} /. μ -> ν, {
DirichletCondition[{u[x, y] == U0[y, i*t0], v[x, y] == 0},
x == 0.],
DirichletCondition[{u[x, y] == 0., v[x, y] == 0.},
0 <= x <= L && y == 0 || y == H],
DirichletCondition[{u[x, y] == 0,
v[x, y] == 0}, (x - x0)^2 + (y - y0)^2 == r0^2],
DirichletCondition[p[x, y] == P0[i - 1][x, y], x == L]}}, {u, v,
p}, {x, y} ∈ Ω,
Method -> {"FiniteElement",
"InterpolationOrder" -> {u -> 2, v -> 2, p -> 1},
"MeshOptions" -> {"MaxCellMeasure" -> 0.001}}], {i, 1, K}];
{ContourPlot[UX[K/2][x, y], {x, y} ∈ Ω,
AspectRatio -> Automatic, ColorFunction -> "BlueGreenYellow",
FrameLabel -> {x, y}, PlotLegends -> Automatic, Contours -> 20,
PlotPoints -> 25, PlotLabel -> u, MaxRecursion -> 2],
ContourPlot[VY[K/2][x, y], {x, y} ∈ Ω,
AspectRatio -> Automatic, ColorFunction -> "BlueGreenYellow",
FrameLabel -> {x, y}, PlotLegends -> Automatic, Contours -> 20,
PlotPoints -> 25, PlotLabel -> v, MaxRecursion -> 2,
PlotRange -> All]} // Quiet
{DensityPlot[UX[K][x, y], {x, y} ∈ Ω,
AspectRatio -> Automatic, ColorFunction -> "BlueGreenYellow",
FrameLabel -> {x, y}, PlotLegends -> Automatic, PlotPoints -> 25,
PlotLabel -> u, MaxRecursion -> 2],
DensityPlot[VY[K][x, y], {x, y} ∈ Ω,
AspectRatio -> Automatic, ColorFunction -> "BlueGreenYellow",
FrameLabel -> {x, y}, PlotLegends -> Automatic, PlotPoints -> 25,
PlotLabel -> v, MaxRecursion -> 2, PlotRange -> All]} // Quiet
dPl = Interpolation[
Table[{i*t0, (P0[i][.15, .2] - P0[i][.25, .2])}, {i, 0, K, 1}]];
cD = Table[{t0*i, NIntegrate[(-ν*(-Sin[θ] (Sin[θ]
\!\(\*SuperscriptBox[\(UX[i]\),
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x0 + r Cos[θ],
y0 + r Sin[θ]] + Cos[θ]
\!\(\*SuperscriptBox[\(UX[i]\),
TagBox[
RowBox[{"(",
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x0 + r Cos[θ],
y0 + r Sin[θ]]) + Cos[θ] (Sin[θ]
\!\(\*SuperscriptBox[\(VY[i]\),
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x0 + r Cos[θ],
y0 + r Sin[θ]] + Cos[θ]
\!\(\*SuperscriptBox[\(VY[i]\),
TagBox[
RowBox[{"(",
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x0 + r Cos[θ],
y0 + r Sin[θ]]))*Sin[θ] -
P0[i][x0 + r Cos[θ], y0 + r Sin[θ]]*
Cos[θ]) /. {r -> r0}, {θ, 0, 2*Pi}]}, {i,
1000, 2000}]; // Quiet
cL = Table[{t0*i, -NIntegrate[(-ν*(-Sin[θ] (Sin[θ]
\!\(\*SuperscriptBox[\(UX[i]\),
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x0 + r Cos[θ],
y0 + r Sin[θ]] + Cos[θ]
\!\(\*SuperscriptBox[\(UX[i]\),
TagBox[
RowBox[{"(",
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x0 + r Cos[θ],
y0 + r Sin[θ]]) +
Cos[θ] (Sin[θ]
\!\(\*SuperscriptBox[\(VY[i]\),
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x0 + r Cos[θ],
y0 + r Sin[θ]] + Cos[θ]
\!\(\*SuperscriptBox[\(VY[i]\),
TagBox[
RowBox[{"(",
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x0 + r Cos[θ],
y0 + r Sin[θ]]))*Cos[θ] +
P0[i][x0 + r Cos[θ], y0 + r Sin[θ]]*
Sin[θ]) /. {r -> r0}, {θ, 0, 2*Pi}]}, {i,
1000, 2000}]; // Quiet
{ListLinePlot[cD,
AxesLabel -> {"t", "\!\(\*SubscriptBox[\(c\), \(D\)]\)"}],
ListLinePlot[cL,
AxesLabel -> {"t", "\!\(\*SubscriptBox[\(c\), \(L\)]\)"}],
Plot[dPl[x], {x, 0, 8}, AxesLabel -> {"t", "ΔP"}]}
f002 = FindFit[cL, a*.5 + b*.8*Sin[k*16*t + c*1.], {a, b, k, c}, t]
Plot[Evaluate[a*.5 + b*.8*Sin[k*16*t + c*1.] /. f002], {t, 4, 8},
Epilog -> Map[Point, cL]]
k0=k/.f002;
Struhalnumber = .1*16*k0/2/Pi
cLm = MaximalBy[cL, Last]
sol = {Max[cD[[All, 2]]], Max[cL[[All, 2]]], Struhalnumber,
dPl[cLm[[1, 1]] + Pi/(16*k0)]}
In Fig. 1 shows the components of the flow velocity and the required coefficients. Our solution of the problem and what is required in the test
{3.17805, 1.03297, 0.266606, 2.60427}
lowerbound= { 3.2200, 0.9900, 0.2950, 2.4600};
upperbound = {3.2400, 1.0100, 0.3050, 2.5000};
![Fig1][2]
Note that our results differ from allowable by several percent, but if you look at all the results of Table 4 from the cited article, then the agreement is quite acceptable.The worst prediction is for the Strouhal number. We note that we use the explicit Euler method, which gives an underestimate of the Strouhal number, as follows from the data in Table 4.
The next test differs from the previous one in that the input speed varies according to the `U0[y_, t_] := 4*Um*y/H*(1 - y/H)*Sin[Pi*t/8]`. It is necessary to determine the time dependence of the drag and lift parameters for a half-period of oscillation, as well as the pressure drop at the last moment of time. In Fig. 2 shows the components of the flow velocity and the required coefficients. Our solution of the problem and what is required in the test
sol = {3.0438934441256595`,
0.5073345082785012`, -0.11152933279750943`};
lowerbound = {2.9300, 0.4700, -0.1150};
upperbound = {2.9700, 0.4900, -0.1050};
![Fig2][3]
For this test, the agreement with the data in Table 5 is good. Consequently, the two tests are almost completely passed.
I wrote and debugged this code using Mathematics 11.01. But when I ran this code using Mathematics 11.3, I got strange pictures, for example, the disk is represented as a hexagon, the size of the area is changed.
![Fig3][4]
In addition, the numerical solution of the problem has changed, for example, test 2D2
{3.17805, 1.03297, 0.266606, 2.60427} v11.01
{3.15711, 1.11377, 0.266043, 2.54356} v11.03
The attached file contains the working code for test 2D3 describing the flow around the cylinder in a flat channel with a change in the flow velocity.
[1]: http://community.wolfram.com//c/portal/getImageAttachment?filename=test2D2.png&userId=1218692
[2]: http://community.wolfram.com//c/portal/getImageAttachment?filename=test2D2.png&userId=1218692
[3]: http://community.wolfram.com//c/portal/getImageAttachment?filename=test2D3.png&userId=1218692
[4]: http://community.wolfram.com//c/portal/getImageAttachment?filename=Math11.3.png&userId=1218692
[331]: http://community.wolfram.com//c/portal/getImageAttachment?filename=CylinderRe100test2D2.gif&userId=1218692
[332]: http://community.wolfram.com//c/portal/getImageAttachment?filename=2D2test.png&userId=1218692Alexander Trounev2018-08-31T11:44:04ZUnit Checking in System Modeler and making use of unit attributes
http://community.wolfram.com/groups/-/m/t/1451968
Unit checking (or dimensional analysis) is an important means of model verification:
[Units of Measurement in a Modelica Compiler][1]
**Question 1: Is unit checking implemented in the current version of Wolfram System Modeler?**
(If not, what options are there to achieve it by other means?)
In Modelica units are stored as attributes to variables (e.g. Real( unit = "m/s" ) ) and are usually enabled by modifications of existing basic types with default values for these attributes (e.g. the default unit for the type Real is " ").
**Question 2: How to make use of this "meta information" within the model itself:**
- How can I access the unit of some variable?
parameter String referencedUnit = <function of another variable to access its unit>
- How can I assign such a referenced unit to another variable within the runtime of a model? (e.g. I have written a converter function to convert units of time, say from minutes to months, as this is frequently needed in economic or business simulation models where SIunits are rather an exception).
output convertedTime =Real( quantity = "Time", unit = <to be assigned in func-body> )
or
input String timeUnitA; // unit of time to convert from
input String timeUnitB; //unit of time to convert to
output convertedTime = Real( quantity = "Time", unit = timeUnitB )
(All the above will not work for me as different degrees of variability prohibit this: a modification of the unit attribute must be a constant...).
EDIT: I have cross-posted the second part of my question on stack overflow [(5363743)][2]
[1]: http://lup.lub.lu.se/luur/download?func=downloadFile&recordOId=8889293&fileOId=8889294 "Units of Measurement in a Modelica Compiler"
[2]: https://stackoverflow.com/q/52312562/5363743Guido Wolf Reichert2018-09-12T10:34:04ZMusic Generation with GAN MidiNet
http://community.wolfram.com/groups/-/m/t/1435251
I generate a music with reference to [MidiNet][1]. Most neural network models for music generation use recurrent neural networks. However, MidiNet use convolutional neural networks.
There are three models in MidiNet. Model 1 is Melody generator, no chord condition. Model 2,3 are Melody generators with chord condition. I try Model 1, because it is most interesting in the three models compared in the paper.
**Get MIDI data**
-----------------------
My favorite Jazz bassist is [Jaco Pastorius][2]. I get MIDI data from [here][3]. For example, I get MIDI data of "The Chicken".
url = "http://www.midiworld.com/download/1366";
notes = Select[Import[url, {"SoundNotes"}], Length[#] > 0 &];
There are some styles in the notes. I get base style from them.
notes[[All, 3, 3]]
Sound[notes[[1]]]
![enter image description here][4]
![enter image description here][5]
I change MIDI data to Image data. I fix the smallest note unit to be the sixteenth note. I divide the MIDI data into the sixteenth note period and select the sound found at the beginning of each period. And the pitch of SoundNote function is from 1 to 128. So, I change one bar to grayscale image(h=128*w=16).
First, I create the rule to change each note pitch(C-1,...,G9) to number(1,...,128), C4 -> 61.
codebase = {"C", "C#", "D", "D#", "E" , "F", "F#", "G", "G#" , "A",
"A#", "B"};
num = ToString /@ Range[-1, 9];
pitch2numberrule =
Take[Thread[
StringJoin /@ Reverse /@ Tuples[{num, codebase}] ->
Range[0, 131] + 1], 128]
![enter image description here][6]
Next, I change each bar to image (h = 128*w = 16).
tempo = 108;
note16 = 60/(4*tempo); (* length(second) of 1the sixteenth note *)
select16[snlist_, t_] :=
Select[snlist, (t <= #[[2, 1]] <= t + note16) || (t <= #[[2, 2]] <=
t + note16) || (#[[2, 1]] < t && #[[2, 2]] > t + note16) &, 1]
selectbar[snlist_, str_] :=
select16[snlist, #] & /@ Most@Range[str, str + note16*16, note16]
selectpitch[x_] := If[x === {}, 0, x[[1, 1]]] /. pitch2numberrule
pixelbar[snlist_, t_] := Module[{bar, x, y},
bar = selectbar[snlist, t];
x = selectpitch /@ bar;
y = Range[16];
Transpose[{x, y}]
]
imagebar[snlist_, t_] := Module[{image},
image = ConstantArray[0, {128, 16}];
Quiet[(image[[129 - #[[1]], #[[2]]]] = 1) & /@ pixelbar[snlist, t]];
Image[image]
]
soundnote2image[soundnotelist_] := Module[{min, max, data2},
{min, max} = MinMax[#[[2]] & /@ soundnotelist // Flatten];
data2 = {#[[1]], #[[2]] - min} & /@ soundnotelist;
Table[imagebar[data2, t], {t, 0, max - min, note16*16}]
]
(images1 = soundnote2image[notes[[1]]])[[;; 16]]
![enter image description here][7]
**Create the training data**
-----------------------
First, I drop images1 to an integer multiple of the batch size. Its length is 128 bars and about 284 seconds with a batch size of 16.
batchsize = 16;
getbatchsizeimages[i_] := i[[;; batchsize*Floor[Length[i]/batchsize]]]
imagesall = Flatten[Join[getbatchsizeimages /@ {images1}]];
{Length[imagesall], Length[imagesall]*note16*16 // N}
![enter image description here][8]
MidiNet proposes a novel conditional mechanism to use music from the previous bar to condition the generation of the present bar to take into account the temporal dependencies across a different bar. So, each training data of MidiNet (Model 1: Melody generator, no chord condition) consists of three "noise", "prev", "Input". "noise" is a 100-dimensions random vector. "prev" is an image data(1*128*16) of the previous bar. "Input" is an image data(1*128*16) of the present bar. The first "prev" of each batch is all 0.
I generate training data with a batch size of 16 as follows.
randomDim = 100;
n = Floor[Length@imagesall/batchsize];
noise = Table[RandomReal[NormalDistribution[0, 1], {randomDim}],
batchsize*n];
input = ArrayReshape[ImageData[#], {1, 128, 16}] & /@
imagesall[[;; batchsize*n]];
prev = Flatten[
Join[Table[{{ConstantArray[0, {1, 128, 16}]},
input[[batchsize*(i - 1) + 1 ;; batchsize*i - 1]]}, {i, 1, n}]],
2];
trainingData =
AssociationThread[{"noise", "prev",
"Input"} -> {#[[1]], #[[2]], #[[3]]}] & /@
Transpose[{noise, prev, input}];
**Create GAN**
-----------------------
I create generator with reference to MidiNet.
generator = NetGraph[{
1024, BatchNormalizationLayer[], Ramp, 256,
BatchNormalizationLayer[], Ramp, ReshapeLayer[{128, 1, 2}],
DeconvolutionLayer[64, {1, 2}, "Stride" -> {2, 2}],
BatchNormalizationLayer[], Ramp,
DeconvolutionLayer[64, {1, 2}, "Stride" -> {2, 2}],
BatchNormalizationLayer[], Ramp,
DeconvolutionLayer[64, {1, 2}, "Stride" -> {2, 2}],
BatchNormalizationLayer[], Ramp,
DeconvolutionLayer[1, {128, 1}, "Stride" -> {2, 1}],
LogisticSigmoid,
ConvolutionLayer[16, {128, 1}, "Stride" -> {2, 1}],
BatchNormalizationLayer[], Ramp,
ConvolutionLayer[16, {1, 2}, "Stride" -> {1, 2}],
BatchNormalizationLayer[], Ramp,
ConvolutionLayer[16, {1, 2}, "Stride" -> {1, 2}],
BatchNormalizationLayer[], Ramp,
ConvolutionLayer[16, {1, 2}, "Stride" -> {1, 2}],
BatchNormalizationLayer[], Ramp, CatenateLayer[],
CatenateLayer[], CatenateLayer[],
CatenateLayer[]}, {NetPort["noise"] ->
1, NetPort["prev"] -> 19,
19 -> 20 ->
21 -> 22 -> 23 -> 24 -> 25 -> 26 -> 27 -> 28 -> 29 -> 30,
1 -> 2 -> 3 -> 4 -> 5 -> 6 -> 7, {7, 30} -> 31,
31 -> 8 -> 9 -> 10, {10, 27} -> 32,
32 -> 11 -> 12 -> 13, {13, 24} -> 33,
33 -> 14 -> 15 -> 16, {16, 21} -> 34, 34 -> 17 -> 18},
"noise" -> {100}, "prev" -> {1, 128, 16}
]
![enter image description here][9]
I create discriminator which does not have BatchNormalizationLayer and LogisticSigmoid, because I use [Wasserstein GAN][10] easy to stabilize the training.
discriminator = NetGraph[{
ConvolutionLayer[64, {89, 4}, "Stride" -> {1, 1}], Ramp,
ConvolutionLayer[64, {1, 4}, "Stride" -> {1, 1}], Ramp,
ConvolutionLayer[16, {1, 4}, "Stride" -> {1, 1}], Ramp,
1},
{1 -> 2 -> 3 -> 4 -> 5 -> 6 -> 7}, "Input" -> {1, 128, 16}
]
![enter image description here][11]
I create Wasserstein GAN network.
ganNet = NetInitialize[NetGraph[<|"gen" -> generator,
"discrimop" -> NetMapOperator[discriminator],
"cat" -> CatenateLayer[],
"reshape" -> ReshapeLayer[{2, 1, 128, 16}],
"flat" -> ReshapeLayer[{2}],
"scale" -> ConstantTimesLayer["Scaling" -> {-1, 1}],
"total" -> SummationLayer[]|>,
{{NetPort["noise"], NetPort["prev"]} -> "gen" -> "cat",
NetPort["Input"] -> "cat",
"cat" ->
"reshape" -> "discrimop" -> "flat" -> "scale" -> "total"},
"Input" -> {1, 128, 16}]]
![enter image description here][12]
**NetTrain**
-----------------------
I train by using the training data created before. I use RMSProp as the method of NetTrain according to Wasserstein GAN. It take about one hour by using GPU.
net = NetTrain[ganNet, trainingData, All, LossFunction -> "Output",
Method -> {"RMSProp", "LearningRate" -> 0.00005,
"WeightClipping" -> {"discrimop" -> 0.01}},
LearningRateMultipliers -> {"scale" -> 0, "gen" -> -0.2},
TargetDevice -> "GPU", BatchSize -> batchsize,
MaxTrainingRounds -> 50000]
![enter image description here][13]
**Create MIDI**
-----------------------
I create image data of 16 bars by using generator of trained network.
bars = {};
newbar = Image[ConstantArray[0, {1, 128, 16}]];
For[i = 1, i < 17, i++,
noise1 = RandomReal[NormalDistribution[0, 1], {randomDim}];
prev1 = {ImageData[newbar]};
newbar =
NetDecoder[{"Image", "Grayscale"}][
NetExtract[net["TrainedNet"], "gen"][<|"noise" -> noise1,
"prev" -> prev1|>]];
AppendTo[bars, newbar]
]
bars
![enter image description here][14]
I select only the pixel having the max value among each column of the image, because there is a feature that the image generated by Wasserstein GAN is blurred. I clear the images.
clearbar[bar_, threshold_] := Module[{i, barx, col, max},
barx = ConstantArray[0, {128, 16}];
col = Transpose[bar // ImageData];
For[i = 1, i < 17, i++,
max = Max[col[[i]]];
If[max >= threshold,
barx[[First@Position[col[[i]], max, 1], i]] = 1]
];
Image[barx]
]
bars2 = clearbar[#, 0.1] & /@ bars
![enter image description here][15]
I change the image to SoundNote. I concatenate the same continuous pitches.
number2pitchrule = Reverse /@ pitch2numberrule;
images2soundnote[img_, start_] :=
SoundNote[(129 - #[[2]]) /.
number2pitchrule, {(#[[1]] - 1)*note16, #[[1]]*note16} + start,
"ElectricBass", SoundVolume -> 1] & /@
Sort@(Reverse /@ Position[(img // ImageData) /. (1 -> 1.), 1.])
snjoinrule = {x___, SoundNote[s_, {t_, u_}, v_, w_],
SoundNote[s_, {u_, z_}, v_, w_], y___} -> {x,
SoundNote[s, {t, z}, v, w], y};
I generate music and attach its mp3 file.
Sound[Flatten@
MapIndexed[(images2soundnote[#1, note16*16*(First[#2] - 1)] //.
snjoinrule) &, bars2]]
![enter image description here][16]
**Conclusion**
-----------------------
I try music generation with GAN. I am not satisfied with the result. I think that the causes are various, poor training data, poor learning time, etc.
Jaco is gone. I hope Neural Networks will be able to express Jaco's base.
[1]: https://arxiv.org/abs/1703.10847
[2]: https://en.wikipedia.org/wiki/Jaco_Pastorius
[3]: http://www.bock-for-pastorius.de/midi.htm
[4]: http://community.wolfram.com//c/portal/getImageAttachment?filename=317901.jpg&userId=1013863
[5]: http://community.wolfram.com//c/portal/getImageAttachment?filename=567502.jpg&userId=1013863
[6]: http://community.wolfram.com//c/portal/getImageAttachment?filename=476803.jpg&userId=1013863
[7]: http://community.wolfram.com//c/portal/getImageAttachment?filename=744004.jpg&userId=1013863
[8]: http://community.wolfram.com//c/portal/getImageAttachment?filename=586405.jpg&userId=1013863
[9]: http://community.wolfram.com//c/portal/getImageAttachment?filename=707106.jpg&userId=1013863
[10]: https://arxiv.org/abs/1701.07875
[11]: http://community.wolfram.com//c/portal/getImageAttachment?filename=435507.jpg&userId=1013863
[12]: http://community.wolfram.com//c/portal/getImageAttachment?filename=170508.jpg&userId=1013863
[13]: http://community.wolfram.com//c/portal/getImageAttachment?filename=324809.jpg&userId=1013863
[14]: http://community.wolfram.com//c/portal/getImageAttachment?filename=965210.jpg&userId=1013863
[15]: http://community.wolfram.com//c/portal/getImageAttachment?filename=706311.jpg&userId=1013863
[16]: http://community.wolfram.com//c/portal/getImageAttachment?filename=177112.jpg&userId=1013863Kotaro Okazaki2018-09-02T02:30:04ZReverse the axes of a plot?
http://community.wolfram.com/groups/-/m/t/1459957
Hello and thanks for your help.
I am trying to invert the axes provided by the Plot [] command, to invert the Y axis (vertical) and the graphical maintenance of the x axis (horizontal). Thank you very much for your help, I tried to find an answer in the program itself but I did not find it.
Thank you very much for any help you can give me.Miguel Saldias2018-09-15T19:20:21ZScatter plotting satellites?
http://community.wolfram.com/groups/-/m/t/1454309
Im dealing with a table of GEO satellites, I have generated a table with Az and EL values relative to my location. To bad they couldn't be found in the SatelliteData[] query...
This is a snippet of the data.
dataSatTable = {
{"SAT NAME", "EL", "AZ"},
{"NSS-806", 3.26, 99.47},
{"Galaxy-17-19", 5.69, 258.52},
{"Eutelsat-113", 10.4, 254.51}
}
I need to plot each satellite as a dot with a text tag on a scatter plot, This will form an arc. I then need to add another table of data containing obstructions on the plot.
Any pointers?Mathison Ott2018-09-13T07:18:01Zpsfrag for Mathematica 10
http://community.wolfram.com/groups/-/m/t/474155
Hello,
as far as I understand, psfrag is no longer working with Mathematica 10. Does anyone have a solution for this problem or
knows whether there will be a solution in the near future? Or is there an alternative?
What I want to do is export eps files from Mathematica and include them into Latex with nice Labels.a b2015-04-05T14:34:56ZAvoid issue while using DMSList function?
http://community.wolfram.com/groups/-/m/t/1461490
DMSList[20.365]
The above code should have returned the list of degree, minutes and seconds. But only the following code returned the required results of {20, 21, 54.}.
DMSList[{20.365,0,0}]htan aungmin2018-09-17T09:35:33ZCreate a "Great Circle" on a globe through two given points?
http://community.wolfram.com/groups/-/m/t/1460856
A collegue of mine is on holiday from Amsterdam to Miami. Just for fun I would like to plot the great circle through the center of the earth going through Amsterdam and Miami.
I tried to do that with GeoGraphics/Geopath, but for both came the Error-Message “GeoGraphics/Geopath” is not a graphics primitive. I tried several things but I cannot get it right ? How to create it ?
The code is from the Wolfram help with some adaption. See att.
Thank youChiel Geeraert2018-09-16T17:35:52ZAvoid issues on demonstrations running on a web browser?
http://community.wolfram.com/groups/-/m/t/1457702
I have noticed at least a couple of demonstrations that one can't run from within a web browser (constrained optimization and union bound probability). When accessing these demonstrations and the cursor is placed in the graphic, a notice pops up that states "this demonstration is optimized for desktop". If the CDF file is downloaded to an iPad running the Wolfram Player, it seems to run correctly. What is it about these demonstrations that prevents them from running inside a web browser?Mike Luntz2018-09-14T13:22:00ZDefine an implicit function?
http://community.wolfram.com/groups/-/m/t/1457290
Hi there,
I have a theoretic decision model that involves a number of equations that cannot be solved in a closed form, i.e. the solution is given only implicitely. The first of these functions then enters another equaition, that again can only be solved implicitely. And encounter problems when I try to implement this. I only need a numerical solution, as it is for illustration prupose only.
Now here's the first equation:
b[x_] := Solve[b == x - 0.1*Quantile[NormalDistribution[0, 1], b], b]
which gives the error message
$RecursionLimit::reclim2: Recursion depth of 1024 exceeded during evaluation of -0.1 Quantile[NormalDistribution,b].
I'd prefer to have the 0.1 as a free parameter "b[x_,s_] :=...", but would be able to live with that if I have to. Any ideas or comments?
Best ChristianChristan Bauer2018-09-14T11:24:04ZConnect Mathematica to a Bluetooth 4.0 device?
http://community.wolfram.com/groups/-/m/t/1454347
How does Mathematica connect to heart rate device via Bluetooth 4.0 and analyze the sampled data? The FindDevices command has tried but the device cannot found. Are there any things to be aware or other methods?
FindDevices[]
{DeviceObject[{"Camera", 1}], DeviceObject[{
"FunctionDemo", 1}], DeviceObject[{
"RandomSignalDemo", 1}], DeviceObject[{"WriteDemo", 1}]}Tsai Ming-Chou2018-09-13T09:37:35ZSolve the following equations with floor and ceil using Wolfram|Alpha?
http://community.wolfram.com/groups/-/m/t/1457369
When I enter this:
n=50, B=150, k=29.5, h=13.5, x=floor(B/k), H=(ceil(n/(2*x))*h)
I get a solution:
> B = 150, h = 27/2, H = 135/2, k = 59/2, n = 50, x = 5
BUT when I enter this:
n=50, B=150, R=2, k=29.5, h=13.5, x=floor(B/k), H=(ceil(n/(R*x))*h)
I get an error:
> Wolfram|Alpha doesn't understand your query
Why?Cev Ing2018-09-14T11:30:28ZIssues with old postings
http://community.wolfram.com/groups/-/m/t/1455753
Are there any issues with old postings in the community? When trying to access an old post, I'm unable to view the page.
Example
Simple Pendulum Experiment Using Mathematica's Image Processing Capability
We can use the pin notebook to capture the position and time in space of the pendulum.[mcode]pics = {}; Pause[10]; Do[AppendTo[pics, {AbsoluteTime[], CurrentImage[]}]; Pause[0.1], {i,...
AUTHOR: Diego Zviovich 5221Views 0replies
HyperLink
http://community.wolfram.com/groups/-/m/t/193779?p_p_auth=xTZ8jYHy
Result:
GROUPS:
Be respectful. Review our Community Guidelines to understand your role and responsibilities. Community Terms of UseDiego Zviovich2018-09-13T23:03:06ZConvert Wolfram Dataset to JSON or CSV via API?
http://community.wolfram.com/groups/-/m/t/1455513
I am getting wolfram CDF format from this,
beta = APIFunction[{"tablename" -> "String"},ResourceData[ResourceObject[#tablename] ]& ]
co = CloudDeploy[beta, Permissions->"Public"]
Response:
Dataset[{<|"Name" -> "Aachen", "ID" -> "1", "NameType" -> "Valid", "Classification" -> "L5", "Mass" -> Quantity[21, "Grams"], "Fall" -> "Fell", "Year" -> DateObject[{1880}, "Year", "Gregorian", -5.], "Coordinates" -> GeoPosition[{50.775, 6.08333}]|>, <|"Name" -> "Aarhus", "ID" -> "2", "NameType" -> "Valid", "Classification" -> "H6", "Mass" -> Quantity[720, "Grams"], "Fall" -> "Fell", "Year" -> DateObject[{1951}, "Year", "Gregorian", -5.], "Coordinates" -> GeoPosition[{56.18333, 10.23333}]|>, <|"Name" -> "Abee", "ID" -> "6", "NameType" -> "Valid", "Classification" -> "EH4", "Mass" -> Quantity[107000, "Grams"], "Fall" -> "Fell", "Year" -> DateObject[{1952}, "Year", "Gregorian", -5.], "Coordinates" -> GeoPosition[{54.21667, -113.}]|>, <|"Name" -> "Acapulco", "ID" -> "10", "NameType" -> "Valid", "Classification" -> "Acapulcoite", "Mass" -> Quantity[1914, "Grams"], "Fall" -> "Fell", "Year" -> DateObject[{1976}, "Year", "Gregorian", -5.], "Coordinates" -> GeoPosition[{16.88333, -99.9}]|>, <|"Name" -> "Achiras", "ID" -> "370", "NameType" -> "Valid", "Classification" -> "L6", "Mass" -> Quantity[780, "Grams"], "Fall" -> "Fell", "Year" -> DateObject[{1902}, "Year", "Gregorian", -5.], "Coordinates" -> GeoPosition[{-33.16667, -64.95}]|> }]
I need this in a JSON format, I tried to convert it using URLexecute it didn't work
Does anyone know any pythonic or wolfram way to convert this into JSON or CSV?Sag Mk2018-09-13T17:28:29ZObtain an inhomogeneous compound Poisson process?
http://community.wolfram.com/groups/-/m/t/1456738
Mathematica has functions of compound Poisson process and inhomogeneous Poisson process. But it does not have a function of the combination of the two. In other words, is it possible to obtain a compound Poisson process with time-varying intensity?
Thanks in advance!Livvy Zhen2018-09-14T03:36:21ZPrevent the reset of the DynamicImage settings?
http://community.wolfram.com/groups/-/m/t/1455313
Dear all,
In the following code how can I prevent the reset of the definition of the dynamicImage, i.e., if a have zoomed the image and then move the slider how can I prevent that the dynamicImage lose the zoom:
img = ExampleData[{"TestImage", "House"}];
image = {img, img*2.0};
Manipulate[
Row[{DynamicImage@img, image[[ind]]}],
{ind, 1, 2, 1}
]
Please note this example is a toy problem.
Thank you,
LuisLuis Mendes2018-09-13T14:41:01Z