@acer  Thanks!

This works quite well. I'll try and implement it into the rest of the code. 

Thanks again!

@acer thanks for the reply! 

I had trialed a few integration methods, however the function wasn't an ArrayInterpolation in these trials, it was a function given out by dsolve. Anyway, in those trials, I found the 2D integration CubaCuhre to be the best. Anyway, I have created a simplified version of what I am trying to accomplish. 

Thanks again. 

Download SimplifiedIntegration.mwSimplifiedIntegration.mw

@Carl Love 

Thanks for the quick reply. I would like to test this out, however I'm using Maple 18.02, which doesn't seem to incorporate the lowess function. :(

@ecterrab 

Perfect. 

I hadn't tried redefining the g_[] metric. In hindsight it seems an obvious solution. 

This is the solution I was chasing.

Thanks

@acer  

I've had a look at this problem more closely. The issues lies with my stucture of V. In the rho range at hand, I get complex exponentials which are not dealt with correctly. Fundamentally this is an issue with the V function I have defined. I will need to correct this and/or work out how to handle what I am expecting in this complex regime. 

All the numerical solving techniques you have provided will then be very useful. Thanks again for your help and hopefully this thread helps others in the future. 

Cheers

@acer Yes, I was busy trying to get the numerical ODE to solve and hadn't payed attention to the actual IVP I had sent. I'm not sure if you're interested in the problem, but its a physics one. 

I have the function V, which is a potential. It contains three parts V_0, VCWc and V_TDD. Because of the complications of V_TDD, I have given it default values of lambda and mu (defined at the start of the sheet). The other two, I have used more accurate coefficients for lambda and mu, namely the solutions of the fsolve argument where I've set rho=250 and 125^2. Func is essentially V_C + VCW (if mu and lambda were undefined variables). Lets call this "Func". I take the first derivative of of Func wrt to rho and at the point rho=246 this derivative should be zero. Similarly the second derivative should be 125^2. These conditions give me values of mu and lambda as a function of kappa. 

The physical problem is that I am writting a modified potential to the Higgs field. Its derivative at the known minimum (electroweak minumum 246 Gev) should be zero and the second derivative should be the higgs mass, 125^2. 

In regards to accuracy, I'm trying to solve this with what ever accuracy I can get. The ODE I'm solving is like a ball rolling down a hill with friction. The potential looks like 

where the blue line is the negative potential. I need to find the initial rho value which gives a result such that it rolls down the potential and lands onto the second bump. I do this via a bisection. Essentially I pick some initial point p(0.001)=190, solve my ODE see if it goes over or under the bump. Then repeat with a higher/lower value until I narrow down to some degree of accuracy where it sits on the second hill for a little bit before falling back to its minimum (in this case around 0). I'm using the Student[NumericalAnalysis]-Roots procedure to do this. 

This means, I need to repetitively solve this IVP and require speed. I spose I dont need a rediculous level of accuracy for each iteration. 

I will have a play with some polynomial fitting functions to approximate this derivative. 

Thanks again!

@acer 

Hi Acer, thank you VERY much for this amazing response. I've learnt a lot about useful methods for numerical differentiion. 

I just have one further, hopefully small problem. 

I had intentionally used the negative derviative (was used for an entirely different reason), however in this problem, it needs to be positive. The solution I'm looking for ranges over a larger range of rho. Typically for Rho from -50 to 250. 

The solutions presented here work very well for rho from 1 to 10. However the symbolic differentiation breaks down between rho = 50 and 150. (this was one of the main reasons I chose to numerically differentiate). 

The solutions you have presented for mapping a polynomial of degree n is very useful, however I am having difficulty mapping a polynomial (ideally to the dVdrho function you have created, or alternatively to the fdiff function). However I am unable to do this, the mapping to the fdiff function seems to take an unpractically long time, do you know of a method to speed this mapping up?. 

Thanks in advance. 

I've attached your previous worksheet with a few comments on the end. 

Download help2.mw

@acer 

Hi acer, thanks for the indepth response. For some reason I cant enter the contents of this maple file, but here is a link. I've added more complexity (probably too much), but my problem still seems trivial, in that I need to get dsolve to run over some black box function. 

Help.mw

Thanks!

@Kitonum @tomleslie 

A simplified version of the function is essentially this guy:

g:=(z)->Int(Re(x^2*ln(1+exp(-sqrt(x^2+(-.23*z^2)*(1/45000))))), x = 0 .. infinity, digits = 10, method = _d01amc)

I'm using this specific integration technique as its the fastest for solving the full form of the equation. 

@Kitonum, I've tried using fsolve. I can plot the function, and see roughly where the roots lie. Running fsolve(deriv) finds roots which are not accurate and in some cases are blatently wrong.

However you've solved my problem. It was using the ' ' - to wrap around my function. I can now use the Roots package and it gives the correct solutions. Perfect. 

Is there page you could link me to, to help me understand the exact meaning of the ' '?

Thanks guys!!

@ecterrab 

Perhaps I am misunderstanding the math. I am trying to define some tensor e[mu,~a] as you have done above in your second line. I would expect the components of e[mu,~a] to be the same as defined? (Not contracted with the metric). 

The scenario I am picturing is illustrated in the following code (I'm using the same example as before as I've already got the maple code :)). 

Download InverseTetradTensors.mw

@ecterrab 

Hi Edgardo, 

Thanks for the tips. I'm not sure if you want me to keep posting bugs, or just wait a couple of months for the package to be ironed out. 

I'll add just one last one. 

Defining Tensors with tetrad indicies - It doesn't like the ordering or contravariant definitions. 

Eg.

e[mu,~a] = <some matrix>; Define(%)

e[mu,~a,matrix] != <some matrix> but

e[a,~mu, matrix] = <some matrix>

Also

e[~a,mu] = <some matrix>; Define(%)

e[~a,mu,matrix] != <some matrix>

e[~mu,a,matrix] = <some matrix>

In these example's i'm using a non-minkowskian metric. 

@ecterrab 

Hi, I've noticed some odd calculations of Tensors with tetrad indices. Here is a simple example. 

Download TensorArray_Tetrads.mw

@ecterrab 

Thanks for the update. I'll keep updating and let you know if I find any issues or bugs. :)

@ecterrab 

Hi Edgardo, 

The physics package you are developing is amazing. I have looked through the null tetrad worksheet you suggested and it seems to contain everything I was looking for. 

I will stop hasseling you, and let you finish the package in peace, however I was wondering if there was a changelog or something I could subscribe to, in order to keep up to date with your progress. 

One final thing, I was unable to use the Tetrad sub-package, there seems to be a bug or perhaps something wrong with my setup. I'm using the linux version of Maple 18 (fully updated). And I have this:

> with(Physics): <- All good

>with(Physics:-Tetrads); - Error, (in with) user level initialization for package `Physics:-Tetrads' failed: Print is not a command in the Physics package

Oh, I should mention the physics package was updated as of today. 

Thanks again.

@ecterrab 

Hi Edgardo, 

I have seen the SumOverRepeateIndicies (I dont use this as I am summing over different ranges of objects). 

If I have a non-minkwoskian metric, the inverse metric g_[~mu,~nu] is obviously different to the metric g_[mu,nu]. 

I have tensors of this form, where adding the tilde makes life much easier as it automatically does the contraction with the metric for me. 

The command 

> add(g_[~mu,2], ~mu=1..2) does the addition g_[1,2] + g_[2,2] not the desired addition g_[~1,2] + g_[~2,2].  

There may not be an easy way of doing this, and I may be able to reformulate my problem all in terms of tensors.