## 150 Reputation

10 years, 122 days

## 602 integrals to go...

@Axel Vogt Thanks! I will try some more tricks with combine(), Re() etc,. Also, the problem with elminating k1, k2 is that this is not possible for all 603 integrals. I've been playing around with exporting maple to MathML, since mathematica can read in those files. Unfortunetly, I can't find a simple command to export a MathML expession to a .mml file. I tried writeto() which saves as a .txt file? And then tried to read in with MMa's Import[%, "mathML"] with no luck.

Edit: I can convert to Mathematica pretty easily now. I wrote a simple script to replace all "sin(x)" with "Sin[x]", etc....

## yup...

@Markiyan Hirnyk yes, thats what I got too. Perhaps I made a typo when converting from Mma to Maple.

## thanks...

@ecterrab awesome! thanks for the quick reply.

## saved...

@ecterrab Thanks so much! You saved me a world of pain in trying to figure this out!

## Seq...

seq is much much faster than for loops. You should try to re-write the code using only seq. If you have a multi-core computer, you could try Threads:-Seq or Grid:-Seq.

## thanks...

@ecterrab Great thanks for the quick reply!

## No way too save memory?...

@Markiyan Hirnyk So in conclusion there is no way to save memory on this since interpolate has been removed for IVP's?

btw this was related to an earlier post http://www.mapleprimes.com/questions/149602-Solving-First-Order-Linear-Coupled-DE-Matrix-Form where I was trying to solve  a set of first order, linear ODE's in matrix form, I would love to be able to solve all 961 equations for more accuracy, and the minimun is aroudn 132, I can od about 170 before I get a memory error.

X:=Vector([seq(x[i](t),i=1..(L+1)^2)]);
eqs:=convert(Equate(lhs(%),rhs(%)),set);
t0:=1: #Initial t
inits:={seq(x[i](t0)=x0[i],i=1..(L+1)^2)};

res:=dsolve(eqs union inits,numeric,parameters=[seq(x0[i],i=1..(L+1)^2)],range=1..15);
res(parameters=[1,2,3,.1,.2,.3,4,5,6]);
plots:-odeplot(res,[Re(x[1](t)),Im(x[1](t))],1..15,refine=1

 With the range option
 When used with the range option, the method computes the solution for the IVP over the specified range, storing the solution information internally, and uses that information to rapidly interpolate the desired solution value for any call to the returned procedure(s).
 Though possible, it is not recommended that the returned procedure be called for points outside the specified range.
 This method can be used in combination with the refine option of odeplot to produce an adaptive plot (that is,odeplot uses the precomputed points to produce the plot when refine is specified.)
 It is not recommended that this method be used for problems in which the solution can become singular, as each step is stored, and many steps may be taken when near a singularity, so memory usage can become a significant issue.
 The storage of the interpolant in use by this method can be disabled by using the interpolation=false option described below. This is recommended for high accuracy solutions where storage of the interpolant (in addition to the discrete solution) requires too much memory. Disabling the interpolant is not generally recommended because the solution values are obtained from an interpolation of the 5 closest points, and does not necessarily provide an interpolant with order 4 error.

Also the matrix DE is a system of 961 (testing on 161 currently)  linear, coupled, complex valued DEs

## Perfect, thanks! But complex numbers?...

that works surprisingly fast! Thanks! But with the parameters option, can I pass complex numbers? i.e.

coef := x0[1] = .1126537-0.2275632e-2*I, x0[2] = 0.7460254e-1+0.6902865e-1*I, x0[3] = .1339367-1.767657*10^(-7)*I, x0[4] = -0.6367112e-1+0.6881254e-1*I, x0[5] = 0.4677631e-2+0.8336431e-1*I, x0[6] = .1101614+.1061479*I, x0[7] = 0.5578725e-1+0.5915829e-2*I, x0[8] = -0.9790946e-1+.1016116*I, x0[9] = -0.3979998e-2-0.7092869e-1*I

c:=convert(%,list);

res:=dsolve(eqs union inits,numeric,parameters=[seq(x0[i],i=1..(L+1)^2)],range=1..15);
res(parameters=[c] );

## Perfect, thanks! But complex numbers?...

that works surprisingly fast! Thanks! But with the parameters option, can I pass complex numbers? i.e.

coef := x0[1] = .1126537-0.2275632e-2*I, x0[2] = 0.7460254e-1+0.6902865e-1*I, x0[3] = .1339367-1.767657*10^(-7)*I, x0[4] = -0.6367112e-1+0.6881254e-1*I, x0[5] = 0.4677631e-2+0.8336431e-1*I, x0[6] = .1101614+.1061479*I, x0[7] = 0.5578725e-1+0.5915829e-2*I, x0[8] = -0.9790946e-1+.1016116*I, x0[9] = -0.3979998e-2-0.7092869e-1*I

c:=convert(%,list);

res:=dsolve(eqs union inits,numeric,parameters=[seq(x0[i],i=1..(L+1)^2)],range=1..15);
res(parameters=[c] );

## Distinguishing y'(x) & y'...

@mehdi_mech  Yes, if we treat y'=y'(x) its circular to find x in terms of y'(x)  [ y'(x(y'(x....))) doesnt make sense]. But if you merely want everything in terms of a new independant variable, y', which happens to be dy/dx in the previous variable system then we can assign x(y') and y(y') that satisfy the above post it amounts to just a change of variables. I don't think this serves and practical reason but we can obtain x(y') and y(y') up to a numerical constant from y' alone.

## Distinguishing y'(x) & y'...

@mehdi_mech  Yes, if we treat y'=y'(x) its circular to find x in terms of y'(x)  [ y'(x(y'(x....))) doesnt make sense]. But if you merely want everything in terms of a new independant variable, y', which happens to be dy/dx in the previous variable system then we can assign x(y') and y(y') that satisfy the above post it amounts to just a change of variables. I don't think this serves and practical reason but we can obtain x(y') and y(y') up to a numerical constant from y' alone.

## Tentative Yes?...

@Markiyan Hirnyk  given that we found y(y') and y''(y') we can find x(y') from  int( 1/y'' dy' ), but I think my code above also treats x as x(y') since

df/dy' = df/dx * dx/dy  = df/dx *1/y', and we can rearrange and integrate for x(y').

We could do more work and find the inverse function without calculus, if we were given a specific function. i.e. y'=sin(x), x=arcsin(y') is the inverse. So I tentatively say yes, it's  possible and already implemented.

## Tentative Yes?...

@Markiyan Hirnyk  given that we found y(y') and y''(y') we can find x(y') from  int( 1/y'' dy' ), but I think my code above also treats x as x(y') since

df/dy' = df/dx * dx/dy  = df/dx *1/y', and we can rearrange and integrate for x(y').

We could do more work and find the inverse function without calculus, if we were given a specific function. i.e. y'=sin(x), x=arcsin(y') is the inverse. So I tentatively say yes, it's  possible and already implemented.

## Proof that Imaginary Part is Zero, not r...

@Axel Vogt I was pointing out another mistake; sometimes Maple is off by a negitive sign like the (3,-3),(3,3) case and other times it is off mixing up real/img. parts like (0,0),(3,3).

I proved that the imaginary part must be zero...not the real part, can you point out where I might have made a mistake?

What other conventions are there? I'm using the standard physics one found in Griffiths, Mathematica, Wikipedia, etc.

I wrote my own SphericalY proc that fixes the issue (I think) and heres the file that compares their results in the (3,-3),(0,0) case. I noticed they only disagree when I have something of the form Y(L,-L,theta,phi).

bug_in_sphericals.mw

 1 2 3 Page 1 of 3
﻿