It looks like Maple ran out of memory trying to solve your equation symbolically. Get a numeric solution instead.

You can plot the solutions with plots:-odeplot. See ?dsolve,numeric and ?plots,odeplot .

Example: Plot t vs. L(t):

plots:-odeplot(Sol, [[t,L(t)]], t= 0..100)

For some reason I can't compile on my computer. (And I'd appeciate some advice on how to correct that.)

Anyway, Acer's evalhf code runs in 2.55 secs on my computer. By parallelizing that and making some small tweaks to the innermost loop (removing the extra trunc), I got that down to 0.815 secs (modulo my eight processors). The parallelization should also work with the compiled version, but I can't test that. (But I left in the code for the compiled version, of course.) Here's the code: Simu_para_evalhf.mw

Like this:

BVP:= {diff(y(x),x$2)-2*y(x)=0,y(0)=1.2,y(1)=0.9}:
Sol_num:= dsolve(BVP, numeric, method=bvp[midrich], output= listprocedure):
Sol_exact:= dsolve(BVP):
Res:= eval(y(x), Sol_num) - unapply(eval(y(x), Sol_exact), x):
plot(Res, 0..1);

The remaining time can be cut in half (!!!) simply by replacing frac(x) with x - trunc(x) since trunc(x) has already been computed. frac is symbolic, and trunc is kernel. It's amazing how much of a difference this makes.

Like most code that generates a lot of random numbers, this one is an obvious candidate for parallelization. It just required a trivial change to your outermost procedure and converting writable globals to locals (which is good coding practice anyway). The benefit that you derive from this will be proportional to the number of processors that you have. I have eight, most modern home computers have four. Windows Resource Monitor showed that I had 100% utilization on all processors, so parallelization can't get any better than that.

Making these two changes, I got the time down to 6 seconds. Here is the modified code: Simulation_Para.mw

There should be no loop for i from 1 to Size in procedure SampleQ. The loop in Qsim is already for i from 1 to Size. So SampleQ should begin

SampleQ:=proc()
    global Q,obs,forv,A,rSum,cSum;
    local number,i,j,k,n,m,broekdel;    
   
     obs:=Matrix(r,c):   #Creating a r times c matrix with zeroes - to hold the observed values
     forv:=Matrix(r,c):   #Creating a r times c matrix with zeroes - to hold the expected values
     A:=Sample(X,1..Total):   

...

This reduces the run time to 51 seconds for me.

Since I think that you're using a less-than-most-recent version of Maple, see ?grelgroup . If you've upgraded to Maple 17 since we've last talked about it, then see ?GroupTheory,Group . In either case, "words" are represented as lists; so instead of a*b*a^(-1), you'd use [a,b,a^(-1)] or [a,b,1/a].

To represent the group that you used in your example:

G:= grelgroup({a,b}, {[a,a], [b,b], [a,b,1/a,1/b]});

(You can use a^(-1) instead of 1/a if you prefer.)

In this representation, a and b are called the generators.

G is a group, but not a permutation group because its elements are not permutations (a is just a letter, obviously; so it's not a permutation). However, every group is isomorphic to some permutation group. G has four elements: [e, a, b, ab] (where e is the identity element). Make that list correspond to [1,2,3,4]. Now we wish to find the permutation that corresponds to a. Multiply each element of the list by a (this group is abelian (commutative), so order doesn't matter) and you get [ae, aa, ab, aab] = [a, e, ab, b], which corresponds to [2,1,4,3], which is [[1,2],[3,4]] in disjcyc notation. So, in a sense of "is", a "is" the permutation [[1,2],[3,4]], but only because of the order that we listed the elements. Likewise, if we multiply the original list by b, we find that b corresponds to the permutation [[1,3], [2,4]]. Those two elements are enough to generate the group. We can now say

GP:= permgroup(4, {a= [[1,2],[3,4]], b=[[1,3],[2,4]]});

Now, it turns out that this is the minimal-degree permutation representation of this particular group, but that was just luck. The technique outlined here will always give a representation of degree the same as the order of the group, but there is usually a representation of lesser degree.

I hope that that increased your understanding and didn't cause more confusion.

Another way:

[seq](lambda - x^2, x= r);

The  way that looks closest to the original is

lambda -~ r^~2;

I believe these elementwise operators have been available since Maple 13.
 

What you call a "permutation matrix" is not a permutation matrix. Quoting directly from the first sentence of the Wikipedia article "Permutation Matrix"

a permutation matrix is a square binary matrix that has exactly one entry 1 in each row and each column and 0s elsewhere.

The randomize() command should only be used once per session (once per restart). If it is executed twice in rapid succession, it can choose the same seed, since the seed is based on the clock, and Maple's clock is very coarse-grained for today's computers. Here's an example showing how long it can take for that clock to tick:

R:= proc() randomize(); rand() end proc:
X:= R():
for k while X=R() do end do:
X:= R():
for k while X=R() do end do:
k;
                             73112

You say that you need a different seed for each trial? Why? You can set the seed by passing an argument to randomize:

randomize(1);
randomize(2);

etc. Each will generate a different series of random numbers. But if you don't do the subsequent randomizes, you're going to get different random numbers anyway.

Disjoint cycle notation does not represent fixed points, so there are no singletons in that notation.

The command mulperms is for multiplying permutations. Your f and g are not permutations; they are groups of permutations.

fgen:= convert([3,2,4,1], disjcyc);
                          [[1, 3, 4]]
ggen:= convert([2,4,3,1], disjcyc);
                          [[1, 2, 4]]
f:= permgroup(4, {fgen}):
g:= permgroup(4, {ggen}):
group:-mulperms(fgen, ggen);
                        [[1, 3], [2, 4]]
f1:= group:-RandElement(f);
                          [[1, 3, 4]]
g1:= group:-RandElement(g);
                          [[1, 4, 2]]
group:-mulperms(f1,g1);
                          [[1, 3, 2]]

The interface commands will only affect displays made after you change the setting; it will not affect what is already on the screen. This is also true if you change the settings from the Tools -> Options menu and select Apply to Session.

If you want to see more of the internal structure of Maple objects, a good command to use is lprint.

You asked:

I applied the code with the procedure above but it doesnt converge. would you check that whether I understand your point?

It's a little bit more complcated than that because you have to extract the solution procedures returned by dsolve(..., output= listprocedure) and pass them on to `evalf/Int`.

But in order to effectively work on the convergence, you need to see the error messages returned by dsolve, which are usually informative. When you add the code to trap those errors and pass on the messages, plus the code for all the options to dsolve and `evalf/Int`, it becomes substantially more complicated. I also included code to display a plot of the integrand for each iteration. This is very informative. I did not incude the plots in this post in order to save space. To get the plots, set infolevel[F]:= 3 before calling fsolve.

So, here is my code. I think that I have the continuation parameter in a good place now. I am not getting any dsolve errors. But, the integration is slow and the fsolve is not converging. Hopefully some experts will take a look at this. 

In the ODEs, I gave the integral substitute and the continuation parameter unusual names so that they would be easier to find. The integral substitute is _J_J_, and the continuation parameter is _C_C_. The continuation parameter can be moved any place in the system (or even be in multiple locations) or boundary conditions as long as it effectively disappears when its value is 1. The current placement seems to prevent dsolve's convergence problem. However, the fsolve is still not converging.

Download BVP_Int_fsolve.mw

You wrote:

 I Must stated that I think the value of rhocu must be in order of 1e6

That's what I would've thought also! But I tried dozens of initial points and it always seems to converge towards zero, without actually converging. I think that you may have what is known as "catastrophic cancellation" in the integral. That's when a small result is produced by adding a large positive and a large negative term. You lose significant digits when that happens.

Hypergeometric functions are probably not what you'd consider "ordinary", but the following is the best that I can do with this. It might be possible to convert this answer to another function. To make any progress with this problem in Maple, it is necessary to make an assumption such as assuming f > 1. An additional assumption t > 1 will clean up the signums in the answer.

dsolve(D(y)(t) = sqrt((f-1)/(t^f-t))) assuming f>1, t>1;

Your 16 equations can be divided into four independent sets of equations with four variables each. The first set is

{a1*a3 = x1111, a1*a4= x1112, a2*a3 = x1211, a2*a4= x1212}.

The other three sets are completely analogous. It is almost obvious that there is no solution with the x's completely arbitrary because we must have x1212*x1111 = x1112*x1211 since both sides equal a1*a2*a3*a4.

Thank you for posting this very interesting type of ODE. I am not sure of the exact definition of "well-posed", but this problem is solvable.

I hope that the code below is self-explanatory. If there's anything that you (or any reader) does not understand, please ask.

Download BVP_fsolve.mw