Subject: again: nonlinear fit... I've posted some days before and i didn't get an answer. I think my question was not detailed enough and therefor I'll try it again. My experimental data have to be fit to this model: lamkonv[x_, {v_, Dz_, L_}] := 1/(2*v^2*x^3)*((Dz*x/Pi)^0.5*(L*Exp[-L^2/(4*Dz*x)] - (L + v*x) *Exp[-(L - 2*v*x)^2/(4*Dz*x)]) + 0.5*(L^2 + Dz*x)* (Erf[L/(2*(Dz*x)^0.5)] - Erf[(L - 2*v*x)/(2*(Dz*x)^0.5)])) where Dz,v,L are the parameters. (It is a 2-dim. convective-diffusion model) I've tried an own code, nonlinearFit, NonlinearRegress and other optimization algorithms (quasiNewton..) ad I couldn't work it out. 1.The fit runs in a physically wrong range or 2.I get imaginary values for the parameters or 3. I get the nrlnum -message or 4. The fit does nothing Especially the parameter Dz behaves very strange. The value of Dz i about 10^3 times smaller than the other parameters. I've read in a posting before, that this may cause problems. Can anybody help me with that? === Subject: Re: Poles and Complex Factoring Use the GaussianIntegers option Factor[x^2 + 2x + 10, GaussianIntegers -> True] and get ((1 - 3 [ImaginaryI]) + x) ((1 + 3 [ImaginaryI]) + x) yehuda >the complex factoring is wrong >>I know how to calculate the residue of a fuction using Mathematica, but how >> >can I >>use Mathematica to calculate the order of a complex pole? >>It would also be nice for Mathematica to tell me if a particular singularity >> >is an >>essential singularity, removable singularity or a pole...but this is >> >not >>necessary; just icing on the cake. >>Also, is there a way to factor polynomials with imaginary roots? >>Something like: >> Factor[ x^2 + 2x + 10 ] = (x - 1 + 4.5 I)(x - 1 - 4.5 I) >>Peter >>-- >>Birthdays are good for you: A federal funded project has recently >> >determined >>that people with the most number of birthdays will live the >> >longest..... === Subject: Re: Poles and Complex Factoring Use the GaussianIntegers option Factor[ x^2 + 2x + 10,GaussianIntegers->True ] and get ((1-3 I)+x) ((1+3 I)+x) yehuda -------Original Message------- === Subject: Poles and Complex Factoring the complex factoring is wrong >I know how to calculate the residue of a fuction using Mathematica, but how can I >use Mathematica to calculate the order of a complex pole? >It would also be nice for Mathematica to tell me if a particular singularity is an >essential singularity, removable singularity or a pole...but this is not >necessary; just icing on the cake. >Also, is there a way to factor polynomials with imaginary roots? >Something like: > Factor[ x^2 + 2x + 10 ] = (x - 1 + 4.5 I)(x - 1 - 4.5 I) >Peter === Subject: Re: MathGroup /: Descriptive headings enough in my previous exposition. I do not want to give the impression I think this is the best way to handle the group, nor that this is the best tagging scheme possibile. >So now they are going to send a posting and about a day later >get this bureaucratic message? >Then they have to wait another day to get their posting up? Nope. Every post without a tag will be automatically tagged with the short 'uncertain' tag [] and forwarded to the NG without any delay, just as it happens today. The poster will get the bureaucratic message only once. If he intends to conform to the group he will add a tag in his next postings (or whenever he feels the problem deserves a tag). If he does not care, nothing happens, apart from having a [] prepended to his object. [*] Repliers can act the same way: those who feel the message can be categorized will add a tag, those who don't will leave it as it is. Perhaps the last sentence should read: Posts without a tag will be automatically tagged with []. Your post has already been forwarded to the group with this tag; you do not have to send another post now. >Categorization of postings may look great from the far view but think of >how it looks from the viewpoint of the new poster who is deeply immeshed >in his particular problem and isn't thinking about nice schemes of >organizing postings? You are perfectly right. And that's the role of the 'uncertain' tag. Giving posters the freedom not to tag their messages, and giving repliers the freedom to tag their replies, leaving the OP's subject unaltered (this will simplify a later search with Google, for example). And then there are only ten cathegories, big enough to allow for an easy categorization. >And anyway, Steve has enough to contend with as it is. I started from the premise that a form of tagging was to be suggested. Obviously this is only a lazy poster's view. The moderator can (and surely will) have a different perception of the difficulties involved in the implementation of such a scheme. After all the inherent anarchy of the net can only let moderators suppress unwanted behaviors. Forcing an attitude is another story altogether. : ) [*]The use of the wildcard [] tag can also be of aid in removing multiple Re's, and the need for a tag can help the moderator in filtering out part of the spam. cheers, Peltio invalid address in reply-to. crafty demunging required to mail me. === Subject: Re: NonlinearFit problem Just an observation in addition to those who already mentioned the problem with the first data point. Does your function have to be exactly of the form f= r^a Exp[-b r]? The reason I asked is that in looking at the fit, it appeared to me that one could obtain a better fit by allowing a zero offset in the first parameter. For example f = (r-c)^a Exp[-b (r-c)]. See the steps below to produce a comparison graph. I don't know your purposes for the fit. But if you are fitting to some measurements, the offset c could be reasonable if there was some uncertainty in establishing the actual zero point for the measurement. Adam Smith In[1]:= << Statistics`NonlinearFit` < Hi All, > Could anyone give me any suggestion for the specified fitting function > f= r^a Exp[-b r]? > My data point was given below, > data={{0, 1.00002}, {2.31507, 26.4522}, {4.32033, 56.8265}, {6.63539, > 59.6674}, {8.64066, 39.5536}, {10.9557, 21.6862}, {12.961, > 10.1456}, {15.276, 4.39652}} > The following way, > NonlinearFit[data,f,r,{a,b}], gives the error message, > FindFit::njnum: > The Jacobian is not a matrix of numbers at (a,b)={1.,1.}. > How should I do this fitting without the problem? > Feng-Yin Chang, > Institute of Physics,NCTU,Taiwan === Subject: equal distribution of last digits base ten in the primes by b-normality The {1,3,7,9} last digits of the primes modulo 10 equal distribution conjecture has never been proved, but I have a b- normal iteration for it.. What that says is that the modulo ten function is equally spaced over the base ten. This is the same argument that Dr. Bailey used to say that the digits of Pi are equally probable over base 16 using his Pi digits formula. Thus if Bailey's proof is acceptable so is this. So with experimental evidence of several million primes and this type of functional evidence/proof it has been pretty well estsablished that the four last digits appear equally. Clear[x,a,digits,f] (* designed covergent sum and b- normal iterator based on the Prime first digits modulo 10*) (* sorted iterative randoms form a devil's staircase like step *) f[n_]=1/((10-Mod[Prime[n],10])*10^n) digits=200 a=Table[N[f[n],digits],{n,1,digits}]; b=N[Apply[Plus,a],digits] x[n_]:=x[n]=Mod[10*x[n-1]+1/(10-Mod[Prime[n],10]),1] x[0]=0 Clear[a,b] a=Table[N[x[n],digits],{n,0,digits}]; ListPlot[a,PlotJoined->True,PlotRange->All] b=Sort[Table[N[x[n],digits],{n,0,digits}]]; ListPlot[b,PlotJoined->True,PlotRange->All] Respectfully, Roger L. Bagula tftn@earthlink.net, 11759Waterhill Road, Lakeside,Ca 92040-2905,tel: 619-5610814 : alternative email: rlbtftn@netscape.net URL : http://home.earthlink.net/~tftn === Subject: Re: fourier ( FFT ) > Hi > My file has form: > x1 f(x1) > x2 f(x2) .... > How can I use Mathematica to obtain frequencies ? > Fourier[] works on only one list! > When i tested with f(x)=Sin[a x] sometimes > i obtained wrong result. I am not sure what you mean by the wrong result. Ssezi === Subject: Re: fourier ( FFT ) Fourier[] assume that your data equaly spaced sampled on the interval [0,2Pi). You can construct an interpolation function and resample the data. Jens > Hi > My file has form: > x1 f(x1) > x2 f(x2) .... > How can I use Mathematica to obtain frequencies ? > Fourier[] works on only one list! > When i tested with f(x)=Sin[a x] sometimes > i obtained wrong result. === Subject: Re: fourier ( FFT ) I assume the followings: 1) you already have a list of pairs such as myData={{x1,f[x1]},{x2,f[x2]},...} 2) x2-x1=x3-x2=..., that is, equal spaced samples. under this all you need to do is to take the last value of each pair and use Fourier Fourier[Last /@ myData] will do the trick. yehuda >My file has form: >x1 f(x1) >x2 f(x2) .... >How can I use Mathematica to obtain frequencies ? >Fourier[] works on only one list! >When i tested with f(x)=Sin[a x] sometimes >i obtained wrong result. === Subject: Re: Two y-axes in one graph > I want to plot two functions in one graph with two y-axes with > different scales in Mathematica 5.0.1. I found > http://support.wolfram.com/mathematica/graphics/2d/twoaxisgraph.html > but the suggested function TwoAxisPlot doesn't exist: > In[1]:= ?TwoAxisPlot > Information::notfound: Symbol TwoAxisPlot not found. > Do I have to load a package to use it? > Harry Harry, the definition of TwoAxisPlot can be found on the webpage you mentioned above. (Don't forget to < schrieb im Newsbeitrag > I want to plot two functions in one graph with two y-axes with > different scales in Mathematica 5.0.1. I found > http://support.wolfram.com/mathematica/graphics/2d/twoaxisgraph.html > but the suggested function TwoAxisPlot doesn't exist: > In[1]:= ?TwoAxisPlot > Information::notfound: Symbol TwoAxisPlot not found. > Do I have to load a package to use it? > Harry > -- > When you are in it up to your ears, keep your mouth shut. === Subject: Re: Two y-axes in one graph Hi Harry, Please notice that the link http://support.wolfram.com/mathematica/graphics/2d/twoaxisgraph.html given in your post actually DEFINE the TwoAxisPlot function. To do it efficiently they load the Graphics`MultiplPlot` package. You cannot find this function in any of the Mathematica packages. Just copy and paste the example given there. you can use it immediately. yehuda > I want to plot two functions in one graph with two y-axes with >different scales in Mathematica 5.0.1. I found > http://support.wolfram.com/mathematica/graphics/2d/twoaxisgraph.html >but the suggested function TwoAxisPlot doesn't exist: > In[1]:= ?TwoAxisPlot > Information::notfound: Symbol TwoAxisPlot not found. >Do I have to load a package to use it? >Harry === Subject: Re: Two y-axes in one graph You can use Abort[ ]. The print in the following example is for demonstration only. lmaa := Module[ {}, Do[ If[i > 3, Print[exit with value , i];Abort[]], {i, 1, 5} ]; Return[{not from inside}] ] lmaa; yehuda > I want to plot two functions in one graph with two y-axes with >different scales in Mathematica 5.0.1. I found > http://support.wolfram.com/mathematica/graphics/2d/twoaxisgraph.html >but the suggested function TwoAxisPlot doesn't exist: > In[1]:= ?TwoAxisPlot > Information::notfound: Symbol TwoAxisPlot not found. >Do I have to load a package to use it? >Harry === Subject: Fibonacci based sum that is b-normal on binary numbers This sum and it's b-normal sequence is due to work of a friend who doesn't like me to use his name here or elsewhere .He came up with two very nice sums using Fibonacci numbers. I used the Binet function in them and got very good agreement. So I tried them in a b-normal. I had to modify the result some to get this result. I get a new sum that appears irrational and an iteration that is b-normal . I think that using the Binet function in this makes it a new sequence sum. I thought that this was a very remarkable result. Clear[x,a,digits,f,fib] (* convergent sum based on Fibonacci sequence to make a binary b-normal iteration *) digits=200 fib[n_Integer?Positive] :=fib[n] = fib[n-1]+fib[n-2] fib[0]=0;fib[1] = fib[2] = 1; sfib=Sum[fib[n]/((n+1)*2^(n+1)),{n,0,digits}] N[sfib,digits] x[n_]:=x[n]=Mod[2*x[n-1]+fib[n-1]/(2*n),1] x[0]=0 a=Table[N[x[n],digits],{n,0,digits}] ListPlot[a,PlotJoined->True,PlotRange->All] b=Sort[Table[N[x[n],digits],{n,0,digits}]]; ListPlot[b,PlotJoined->True,PlotRange->All] Respectfully, Roger L. Bagula tftn@earthlink.net, 11759Waterhill Road, Lakeside,Ca 92040-2905,tel: 619-5610814 : alternative email: rlbtftn@netscape.net URL : http://home.earthlink.net/~tftn === Subject: List plotting a vector field containing invalid data I'm solving a 2D CFD problem, which is generating a lot of co-located 2D velocity field data, using a regular rectangular grid. Some of the grid is masked out, representing ocean bottom topography. Is there a relatively simple way to modify ListPlotVectorField to ignore the masked out values (which are stored as as unphysically large values, but this is easily modified) === Subject: Akiyama-Thurston [1,5] minimal Pisot tile This tile was proposed by Dr. Thurston in one of his lectures, but only actually found by Dr. Akiyama. The IFS form was found by S. R. Hinsley. (*Akiyama-Thurston [1,5] minimal Pisot tile*) NSolve[x3-x-1==0,x] c=-0.662358978622373051 s=-0.562279512062301289 c1=0.460202188254280750782172413 s1=0.1825822545574430797121307804 (* Wellin IFS program type*) (* Akiyama_Hinsley*) f1[{x_,y_}] = {x*c-s*y+c, s*x+c*y+s}; f2[{x_,y_}] = {x*c1-s1*y+c1, s1*x+c1*y+s1}; f[x_] := Which[(r=Random[]) <= 3/4, f1[x], r <= 1.00, f2[x]] a=NestList[f, {0,0}, 10000]; Respectfully, Roger L. Bagula tftn@earthlink.net, 11759Waterhill Road, Lakeside,Ca 92040-2905,tel: 619-5610814 : alternative email: rlbtftn@netscape.net URL : http://home.earthlink.net/~tftn === Subject: minimal Pisot tile {2,3] type definition in Mathematic There is an IFS rotation procedure used with twin dragon tile to make Heighway dragon tile called a Riddle rotation. I found this tile using that type of experiment on my IFS. Hinsley has the same tile as well. The [1,5] and [2,3] are the powers of the base Minimal Pisot complex number polynomial solution. (* minimal Pisot tile {2,3] type definition in Mathematica*) c=0.868837 a=220.328 r0=c w0=Pi*a/180 x0=r0*Cos[w0] y0=r0*Sin[w0] x5=r03*Cos[3*w0] y5=r03*Sin[3*w0] x3=r02*Cos[2*w0] y3=r02*Sin[2*w0] t=1 aa=(x*x5-y*y5) bb=(x*y5+y*x5) cc=Cos[t*Pi] ss=Sin[t*Pi] x1=aa*cc-bb*ss+x5+(x5)*t y1=aa*ss+bb*cc+y5-(x5)*t (* Wellin IFS program type*) (* Akiyama_23: curley tile*) f1[{x_,y_}] = {x*x3-y*y3+x3, x3*y+y3*x+y3}; f2[{x_,y_}] = {x1, y1}; f[x_] := Which[(r=Random[]) <= 1/2, f1[x], r <= 1.00, f2[x]] ifs[n_] := Show[Graphics[{PointSize[.001], Map[Point, NestList[f, {0,0}, n]]}], PlotRange->All,AspectRatio->Automatic] Respectfully, Roger L. Bagula tftn@earthlink.net, 11759Waterhill Road, Lakeside,Ca 92040-2905,tel: 619-5610814 : alternative email: rlbtftn@netscape.net URL : http://home.earthlink.net/~tftn ------------------------------------------------------------------------ === Subject: Re: First question > How to exit directly and immediately a module from inside a loop with > a certain value ? > Return seems not to work. > Example: > lmaa := Module[ > {}, > Do[ > If[i > 3, Return[{from inside}]], > {i, 1, 5} > ]; > Return[{not from inside}] > lmaa > {not from inside} > Best > Gilbert Gosseyn Read in the docs about Catch[] and Throw[]: lmaa := Module[{}, Do[If[i > 3, Throw[{from inside}]], {i, 1, 5}]; {not from inside}] Catch[lmaa] {from inside} -- Peter Pein Berlin === Subject: Re: First question the manual say Unless an explicit Return is used, the value returned by Do is Null. and that mean tha Return[] in a Do[]-Loop terminate the Do[] loop and not the function. You mean lmaa [] := Module[{i}, res = Do[ If[i > 3, Return[{from inside}]], {i, 0, 5} ]; ] or lmaa [] := Module[{i = 0}, While[i < 6, If[i > 3, Return[{from inside]; i++ ]; {not from inside} ] lmaa [] := Module[{i}, For[i = 0, i < 6, i++, If[i > 3, Return[{from inside}]] ]; {not from inside} ] Jens > How to exit directly and immediately a module from inside a loop with > a certain value ? > Return seems not to work. > Example: > lmaa := Module[ > {}, > Do[ > If[i > 3, Return[{from inside}]], > {i, 1, 5} > ]; > Return[{not from inside}] > lmaa > {not from inside} > Best > Gilbert Gosseyn === Subject: Re: First question Hi Gilbert, > How to exit directly and immediately a module from inside a loop with > a certain value ? > Return seems not to work. > Example: > lmaa := Module[ > {}, > Do[ > If[i > 3, Return[{from inside}]], > {i, 1, 5} > ]; > Return[{not from inside}] You have a bug in your code, as well as some possible misunderstanding, and certainly some dangerous programming. The bug is the If statement which should read If[i > 3, Return[{from inside}], Return[{not from inside}] ] Your final statement means that the function will always return not from inside. There is a confusion since you may have another i which is global that you want to work from instead. The dangerous programming is that i is a local variable to the loop. It always takes values 1, 2, 3, 4, 5. Which in turn means that it returns the last value (there are no Break[]'s in your code) which is 5. You have some reading to do .... Dave. === Subject: Re: First question lmaa := Catch[Do[If[i > 3, Throw@from inside], {i, 1, 5}]; not from inside] lmaa from inside Bobby > How to exit directly and immediately a module from inside a loop with > a certain value ? > Return seems not to work. > Example: > lmaa := Module[ > {}, > Do[ > If[i > 3, Return[{from inside}]], > {i, 1, 5} > ]; > Return[{not from inside}] > lmaa > {not from inside} > Best > Gilbert Gosseyn -- DrBob@bigfoot.com www.eclecticdreams.net === Subject: Re: First question > How to exit directly and immediately a module from inside a loop with > a certain value ? > Return seems not to work. > Example: > lmaa := Module[ > {}, > Do[ > If[i > 3, Return[{from inside}]], > {i, 1, 5} > ]; > Return[{not from inside}] > lmaa > {not from inside} > Best > Gilbert Gosseyn Method 1 (not really recommended) In[70]:= lmaa := Module[ {}, Do[ If[i > 3, Return[Return[{from inside}]]], {i, 1, 5} ]; Return[{not from inside}] ] lmaa {from inside} Method 2 (recommended): lmaa := Module[ {}, Catch[ Do[ If[i > 3, Throw[{from inside}]], {i, 1, 5} ]; Return[{not from inside}]] ] lmaa {from inside} Finally, there is absolutely no point using a Module whiteout any local variables. Andrzej Kozlowski Chiba, Japan http://www.akikoz.net/~andrzej/ http://www.mimuw.edu.pl/~akoz/ === Subject: Re: Re: newbie question DSolve That use of DSolve returns unevaluated in version 5.0.1, with a different error message: DSolve::litarg To avoid possible ambiguity, the arguments of the dependent variable in `1` should literally match the independent variables. Bobby > If you replace DiracDelta[d - x]*y[x] by the equivalent DiracDelta[d - > x]*y[d] then your equation can be solved as follows (with just one minor > error message appearing twice, which can be ignored) > In[1]:= > s = DSolve[{-y[x] + Derivative[2][y][x] == DiracDelta[d - x]*y[d], y[0] > == 0, y[L] == 0}, y[x], x] > DSolve::nvld : The description of the equations appears to be > ambiguous or > invalid. > Out[1]= > {{y[x] -> 1/2*E^(-d - x)*(-((E^(2*x)*((-1 + E^(2*d))*UnitStep[-d] - > (E^(2*d) - E^(2*L))*UnitStep[-d + L])*y[d])/ > (-1 + E^(2*L))) + ((E^(2*L)*(-1 + E^(2*d))*UnitStep[-d] - > (E^(2*d) - E^(2*L))*UnitStep[-d + L])*y[d])/ > (-1 + E^(2*L)) - E^(2*d)*UnitStep[-d + x]*y[d] + > E^(2*x)*UnitStep[-d + x]*y[d])}} > DSolve::nvld : The description of the equations appears to be > ambiguous or > invalid. > Extracting the solution to u[x] > In[4]:= > u[x_] = y[x] /. s[[1]] > you can Plot it, after assigning numeric values to all relevant quantities: > In[6]:= > L = 1; d = 0.5; y[d] = 1; > Plot[u[x], {x, -1, 4}, PlotRange -> {{-1, 5}, {-2, 1}}]; > Hope this hepls > Wolfgang >> Hello all >> I am trying to use DSolve to solve a ode with discontinuity in it (wave >> equation with a viscous damper injected at a location d) >> This is what i am using >> DSolve[{y''[x]-lamda^2*y[x]==DiracDelta[x-d]*y[x],y[0]==0,y[L]== >> =0},y[x],x] >> the problem I am facing is that >> y[x] on the right hand side (next the delta function) varies w.r.t to >> the location >> y[x]==y[x]&& 0<=x<=d >> y[x]==y[L-x]&&d<=x<=L >> I can solve the above equation without the y[x] coupled to the delta >> function >> Pratik Desai >> ps: This is my third attempt at posting my query, I hope this time it >> makes it to the list :) -- DrBob@bigfoot.com www.eclecticdreams.net === Subject: Re: Re: newbie question DSolve (revisited) I tried your suggestion unfortunately, Mathematica gives me another error it is as follows In[10]:= DSolve[{y''[x] - lamda^2*y[x] == DiracDelta[x - d]*y[d], y[0] == 0, y[l] == 0}, y[x], x] (DSolve::litarg), To avoid possible ambiguity, the arguments of the dependent variable in (the equation) should literally match the independent variables. Pratik Desai ----- Original Message ----- === Subject: Re: newbie question DSolve > If you replace DiracDelta[d - x]*y[x] by the equivalent DiracDelta[d - > x]*y[d] then your equation can be solved as follows (with just one minor > error message appearing twice, which can be ignored) > In[1]:= > s = DSolve[{-y[x] + Derivative[2][y][x] == DiracDelta[d - x]*y[d], y[0] > == 0, y[L] == 0}, y[x], x] > DSolve::nvld : The description of the equations appears to be > ambiguous or > invalid. > Out[1]= > {{y[x] -> 1/2*E^(-d - x)*(-((E^(2*x)*((-1 + E^(2*d))*UnitStep[-d] - > (E^(2*d) - E^(2*L))*UnitStep[-d + L])*y[d])/ > (-1 + E^(2*L))) + ((E^(2*L)*(-1 + E^(2*d))*UnitStep[-d] - > (E^(2*d) - E^(2*L))*UnitStep[-d + L])*y[d])/ > (-1 + E^(2*L)) - E^(2*d)*UnitStep[-d + x]*y[d] + > E^(2*x)*UnitStep[-d + x]*y[d])}} > DSolve::nvld : The description of the equations appears to be > ambiguous or > invalid. > Extracting the solution to u[x] > In[4]:= > u[x_] = y[x] /. s[[1]] > you can Plot it, after assigning numeric values to all relevant > quantities: > In[6]:= > L = 1; d = 0.5; y[d] = 1; > Plot[u[x], {x, -1, 4}, PlotRange -> {{-1, 5}, {-2, 1}}]; > Hope this hepls > Wolfgang >> Hello all >> I am trying to use DSolve to solve a ode with discontinuity in it (wave >> equation with a viscous damper injected at a location d) >> This is what i am using >> DSolve[{y''[x]-lamda^2*y[x]==DiracDelta[x-d]*y[x],y[0]==0,y[L]== >> =0},y[x],x] >> the problem I am facing is that >> y[x] on the right hand side (next the delta function) varies w.r.t to >> the location >> y[x]==y[x]&& 0<=x<=d >> y[x]==y[L-x]&&d<=x<=L >> I can solve the above equation without the y[x] coupled to the delta >> function >> Pratik Desai >> ps: This is my third attempt at posting my query, I hope this time it >> makes it to the list :) === Subject: Re: First question I just copied your post and forgot to change the second Return also to Throw. Andrzej > Throw ? Catch ? , are we talking about Baseball ? :-) > Why does the Catch include 'Return[{not from inside}]' ? > Best, > GG > -- snip >> Method 2 (recommended): >> lmaa := Module[ >> {}, >> Catch[ Do[ >> If[i > 3, Throw[{from inside}]], >> {i, 1, 5} >> ]; >> Return[{not from inside}]] >> ] >> lmaa >> {from inside} >> Finally, there is absolutely no point using a Module whiteout any >> local variables. >> Andrzej Kozlowski >> Chiba, Japan >> http://www.akikoz.net/~andrzej/ >> http://www.mimuw.edu.pl/~akoz/ === Subject: Solve and Reduce Hello list, I want to find q as a function of c, q(c), given the following equation: 2500*c^2 - 25*c^3 + 3500*c*q - 320*c^2*q - 1104*c*q^2 - 1152*q^3 == 0 However, each of the following three methods gives different results. I check the Mathematica Book but still cannot figure out why there are lot! (1) Use Reduce In[5]:= q1[c_] = Reduce[{2500*c^2 - 25*c^3 +3500*c*q - 320*c^2*q -1104*c*q^2 -1152*q^3 == 0, c > 0,q > 0}, q] Out[5]= 0> The main problem seems to be that the differential >> equation does not get evaluated in each step but only once. model (as written) has the same value every time, so why do you want it evaluated more than once? I think you wanted it to depend on a, but you've written it with a->0.05 hardwired in. I think this does what you want: f[t_] := Exp[-0.1*t] + Random[Real, {0, 0.1}]; data = {#, f@#} & /@ Table[t, {t, 0, 10, 1}]; {xData, yData} = Transpose@data; Clear[model, square] model[(a_)?NumericQ] := y /. First[NDSolve[ {Derivative[1][y][t] == (-a)*y[t], y[0] == 1}, y, {t, 0, 10}]] square[a_] := (#1 . #1 & )[ yData - model[a] /@ xData] NMinimize[{square[a], a >= 0}, a] {0.00638361, {a -> 0.0818508}} Here's a more direct (but fully equivalent) approach: Clear[f, model, square] model[a_][t_] := Exp[(-a)*t] f[t_] := model[0.1][t] + Random[Real, {0, 0.1}] data = Table[{t, f[t]}, {t, 0, 10, 1}]; {xData, yData} = Transpose[data]; square[a_] := (#1 . #1 & )[ yData - model[a] /@ xData] NMinimize[{square[a], a >= 0}, a] This always underestimates a, however, because the noise you've added is always positive, and the error function doesn't take that into account. You can do better this way: Clear[f, model, square] model[a_][t_] := Exp[-a*t] f[t_] := model[0.1]@t + 0.1Random[] {xData, yData} = Transpose@Table[{t, f@t}, {t, 0, 10, 1}]; square[a_] := #.# &[yData - a/2 - model[a] /@ xData] NMinimize[{square@a, a >= 0}, a] {0.0139239, {a -> 0.0999625}} and this does pretty well with no minimization at all: Clear[f, model, inverse] model[a_, t_] = Exp[-a*t]; inverse[a_]@{x_, y_} := -Log[y - .5a]/x f[t_] := model[0.1, t] + 0.1Random[] data = Table[{t, f@t}, {t, 1, 10}]; Mean[inverse[0.1] /@ data] 0.100146 I didn't put a lot of science into that. I'm just subtracting the average bias and solving for a in y - .5 a == Exp[-a x]. The LHS of that equation is an unbiased estimate of y (given a and x), but the result isn't necessarily an unbiased estimate of a (given y and x). Bobby > I made a little example to illustrate the problem: > I want to fit a model which is stated as a system of differenttial > equations to some data. > (* curve to be fitted *) > f[t_]:=Exp[-0.1*t]+Random[Real,{0,0.1}]; > data={#,Evaluate[f[#]]}&/@Table[t,{t,0,10,1}]; > ListPlot[data,PlotStyle[Rule]{RGBColor[0,0.2,0.8],PointSize[0.02]}]; > (* the parameter 0.1 is to be found *) > (* the model as a differential equations *) > model:=NDSolve[{y'[t]==-a*y[t],y[0]==1}/.a->0.05,y,{t,0,10}]; > (* minimizing sum of squares *) > NMinimize[{squares=Plus@@Table[(Table[{#,y[t]/.Evaluate[model][[1]]/.t->#}&/@ > Table[t,{t, 0, 10, 1}]][[i, 2]]-data[[i,2]])^2, > {i, Length[data]}], a > 0}, a, > Method -> {DifferentialEvolution}, > StepMonitor :> Print[squares, a], > WorkingPrecision -> 8] > The main problem seems to be that the differential equation does not get > evaluated in each step but only once. > So if anybody could make this example work that would be great. > best, > joerg -- DrBob@bigfoot.com www.eclecticdreams.net === Subject: Re: nonlinear programming with differential-algebraic constraints (2) > I made a little example to illustrate the problem: > I want to fit a model which is stated as a system of differenttial > equations to some data. > (* curve to be fitted *) > f[t_]:=Exp[-0.1*t]+Random[Real,{0,0.1}]; > data={#,Evaluate[f[#]]}&/@Table[t,{t,0,10,1}]; > ListPlot[data,PlotStyle[Rule]{RGBColor[0,0.2,0.8],PointSize[0.02]}]; > (* the parameter 0.1 is to be found *) > (* the model as a differential equations *) > model:=NDSolve[{y'[t]==-a*y[t],y[0]==1}/.a->0.05,y,{t,0,10}]; > (* minimizing sum of squares *) > NMinimize[{squares=Plus@@Table[(Table[{#,y[t]/.Evaluate[model][[1]]/.t->#}&/@ > Table[t,{t, 0, 10, 1}]][[i, 2]]-data[[i,2]])^2, > {i, Length[data]}], a > 0}, a, > Method -> {DifferentialEvolution}, > StepMonitor :> Print[squares, a], > WorkingPrecision -> 8] > The main problem seems to be that the differential equation does not get > evaluated in each step but only once. > So if anybody could make this example work that would be great. > best, > joerg Hi J.9arg, I changed the code to let me play more easily with different ranges and counts of points. I changed f to produce a relative error and finally I used FindMinimum instead of NMinimize because I've got no access to version 5 of Mathematica: In[1]:= (*curve to be fitted*) n = 50; f[t_] := Exp[-0.1*t]*Random[Real, {0.9, 1.1}]; data = {#, Evaluate[f[#]]} & /@ (10 Range[0, n]/n); ListPlot[data, PlotStyle -> {RGBColor[0, 0.2, 0.8], PointSize[0.02]}, Axes -> False, Frame -> True]; (*Plot ommitted*) In[5]:= (*the parameter 0.1 is to be found*) (*the model as a differential equations*) (*evaluate only for numeric a, and increment a counter*) model[a_?NumericQ] := (solvecount++; NDSolve[{y'[t] == -a*y[t], y[0] == 1}, y, {t, Sequence @@ First /@ data[[{1, -1}]]}]); In[6]:= solvecount = 0; FindMinimum[(#1 . #1 & )[ (y[#1] /. Evaluate[model[a]][[1]] & ) /@ First /@ data - Last /@ data ], {a, 0, 1}] solvecount Out[7]= {0.0672762, {a -> 0.101174}} Out[8]= 816 As you can see, model[] has been called 816 times. I hope this helped, Peter -- Peter Pein Berlin === Subject: Re: nonlinear programming with differential-algebraic constraints (2) f[t_] := Exp[-0.1*t] + Random[Real, {0, 0.1}]; data = {#, Evaluate[f[#]]} & /@ Table[t, {t, 0, 10, 1}]; model[a_] := model[a]= y /. NDSolve[{y'[t] == -a*y[t], y[0] == 1} /. a -> 0.05, y, {t, 0, 10}][[1]]; fun[a_?NumericQ] := Module[{ff}, ff = model[a]; Plus @@ ((ff[#[[1]]] - #[[2]])^2 & /@ data) ] NMinimize[{fun[a], a > 0}, a, Method -> {DifferentialEvolution}, StepMonitor :> Print[squares, a], WorkingPrecision -> 8] Jens Joerg Schaber schrieb im Newsbeitrag >I made a little example to illustrate the problem: > I want to fit a model which is stated as a system of differenttial > equations to some data. > (* curve to be fitted *) > f[t_]:=Exp[-0.1*t]+Random[Real,{0,0.1}]; > data={#,Evaluate[f[#]]}&/@Table[t,{t,0,10,1}]; > ListPlot[data,PlotStyle[Rule]{RGBColor[0,0.2,0.8],PointSize[0.02]}]; > (* the parameter 0.1 is to be found *) > (* the model as a differential equations *) > model:=NDSolve[{y'[t]==-a*y[t],y[0]==1}/.a->0.05,y,{t,0,10}]; > (* minimizing sum of squares *) > NMinimize[{squares=Plus@@Table[(Table[{#,y[t]/.Evaluate[model][[1]]/.t->#}&/@ > Table[t,{t, 0, 10, 1}]][[i, 2]]-data[[i,2]])^2, > {i, Length[data]}], a > 0}, a, > Method -> {DifferentialEvolution}, > StepMonitor :> Print[squares, a], > WorkingPrecision -> 8] > The main problem seems to be that the differential equation does not get > evaluated in each step but only once. > So if anybody could make this example work that would be great. > best, > joerg === Subject: Re: Poles and Complex Factoring Factor[x^2+2x+10, GaussianIntegers -> True] (x + (1 - 3*I))*(x + (1 + 3*I)) Bob Hanlon === > Subject: Poles and Complex Factoring > the complex factoring is wrong >I know how to calculate the residue of a fuction using Mathematica, but how > can I >use Mathematica to calculate the order of a complex pole? >It would also be nice for Mathematica to tell me if a particular singularity > is an >essential singularity, removable singularity or a pole...but this is > not >necessary; just icing on the cake. >Also, is there a way to factor polynomials with imaginary roots? >Something like: > Factor[ x^2 + 2x + 10 ] = (x - 1 + 4.5 I)(x - 1 - 4.5 I) >Peter === Subject: Re: First question >How to exit directly and immediately a module from inside a loop >with a certain value ? Return seems not to work. >Example: >lmaa := Module[ >{}, > Do[ > If[i > 3, Return[{from inside}]], > {i, 1, 5} > ]; > Return[{not from inside}] >lmaa >{not from inside} Your problem isn't that Return isn't working but that it only causes an exit from one level of control structure. What happens is as follows when i = 4: The Do loop is exited then the next line of code Return[{not from inside}] is executed causing the module to be exited with the results of this line of code. Change your code as follows: maa[m_] := Module[{n = 0}, Do[If[i > m, n = i; Return[{ from inside}]], {i, 1, 5}]; Return[{StringJoin[n is , ToString[n]],not from inside}]] and see maa[3] {n is 4, not from inside} demonstrating Return does work as advertised. If you want to exit from the Do loop and the Module when i becomes 4, you need to test the value of i immediately after exiting the Do loop and exit the Module based on that test. -- === Subject: Re: First question Actually, I was again not paying attention; you do not need either Return or Throw or anything in the second Return position: lmaa := Module[ {}, Catch[ Do[ If[i > 3, Throw[{from inside}]], {i, 1, 5} ]; not from inside] ] is all you need in this kind of code. You do not need Return or anything else to return output in Mathematica, only to exit certain loops. Andrzej Kozlowski > I just copied your post and forgot to change the second Return also to > Throw. > Andrzej >> *This message was transferred with a trial version of CommuniGate(tm) >> Pro* >> Throw ? Catch ? , are we talking about Baseball ? :-) >> Why does the Catch include 'Return[{not from inside}]' ? >> Best, >> GG >> -- snip > Method 2 (recommended): > lmaa := Module[ > {}, > Catch[ Do[ > If[i > 3, Throw[{from inside}]], > {i, 1, 5} > ]; > Return[{not from inside}]] > ] > lmaa > {from inside} > Finally, there is absolutely no point using a Module whiteout any > local variables. > Andrzej Kozlowski > Chiba, Japan > http://www.akikoz.net/~andrzej/ > http://www.mimuw.edu.pl/~akoz/ === Subject: matlab code to mma? [Please send comments on this directly to Sean - moderator] i have received some codes written in matlab, and i would like to convert it to Mathematica. anyone up for the challenge *and* have the time and energy? It will be nice to see line by line explanation of the conversion. following is the code which I'm having problems understanding. sean ------- function result= gillespie(tf, v0) % function result= gillespie(tf, v0) % % runs the Gillespie algorithm for constitutive expression up to % a time tf. the value of v0 may be specified but has a default value % of 0.01 if nargin == 1 v0= 0.01; end; % initial amounts D= 1; M= 0; N= 0; t= 0; % initialize other variables j= 2; nr= 4; % number of reactions % create some large result matrix (this is not essential but % speeds up the program) result= zeros(100000, 4); % store data for first time point result(1,:)= [t D M N]; while t < tf % conversion s(1)= D; s(2)= M; s(3)= N; % store data result(j,:)= [t D M N]; j= j+1; % calculate propensities a= calpropensities(s, v0); a0= sum(a); % generate two random numbers r= rand(1,2); % calculate time of next reaction dt= log(1/r(1))/a0; % calculate which reaction is next for mu= 1:nr if sum(a(1:mu)) > r(2)*a0 break; end; end; % carry out reaction t= t+dt; if mu == 1 M= M+1; elseif mu == 2 M= M-1; elseif mu == 3 N= N+1; elseif mu == 4 N= N-1; end; end; % store data for last time point result(j,:)= [t D M N]; % remove any excess zeros from result result(j+1:end,:)= []; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function a= calpropensities(s, v0) % rates v1= 0.04; d0= log(2)/180; d1= log(2)/3600; % conversion D= s(1); M= s(2); N= s(3); % propensities a(1)= v0*D; a(2)= d0*M; a(3)= v1*M; a(4)= d1*N; __________________________________ Do you Yahoo!? Check out the new Yahoo! Front Page. www.yahoo.com === Subject: Re: Two y-axes in one graph Hi Harry, > but the suggested function TwoAxisPlot doesn't exist: So create it .... It's defined right there on the page. In the event that you see a picture placeholder, copy the following (* from here *) Needs[Graphics`MultipleListPlot`] TwoAxisListPlot[f_List, g_List, frange_, grange_, opts___?OptionQ] := Module[{old, new, scale, fm, fM, gm, gM, newg}, {fm, fM} = frange; {gm, gM} = grange; scale[var_] = ((var - gm)*(fM - fm))/(gM - gm) + fm ; old = AbsoluteOptions[ListPlot[g, Frame -> True, PlotRange -> grange, DisplayFunction -> Identity], FrameTicks][[1, 2, 2]]; new = (Prepend[Rest[#1], scale[First[#1]]] & ) /@ old; newg = Transpose[{Transpose[g][[1]], Map[scale, Transpose[g][[2]], {1, 2}]}]; MultipleListPlot[f, newg, Frame -> True, FrameTicks -> {Automatic, Automatic, None, new}, SymbolStyle -> {{RGBColor[1, 0, 0]}, {RGBColor[0, 0, 1]}}, FrameStyle -> {{}, {RGBColor[1, 0, 0]}, {}, {RGBColor[0, 0, 1]}}, PlotRange -> frange + .05 (fM - fm), opts]] (* to here *) > In[1]:= ?TwoAxisPlot > Information::notfound: Symbol TwoAxisPlot not found. > Do I have to load a package to use it? You have to create it. Dave.