mehdibgh

240 Reputation

6 Badges

7 years, 276 days

MaplePrimes Activity


These are replies submitted by mehdibgh

@tomleslie Thans,The problem was the way I saved the matrix in matlab.

Here still a problem, since the matrix dimension is 5000*5000 it take five minutes for Maple to import. Any suggestion to speed up the transfer time?

Lets simplify the problem more to find a way to solve the problem.

First drop w and pi for simplicity. Since the function H(w) acts as a continuous Heviside function to change the integral boundaries ( within only positive region of w), the other way of computing the integral is to intersect H with w first, and then try to integrate L inside the upper region of the created area between the intersection of two functions H and w in -1<x,y<1.

NULL

restart

JJ := 5:

II := 5:

with(ArrayTools):

W := RandomArray(II, JJ, M):

``

w := add(add(W[i, j]*LegendreP(i-1, x)*LegendreP(j-1, y), i = 1 .. II), j = 1 .. JJ):

f := add(add((LegendreP(i-1, x)*LegendreP(j-1, y))^2, i = 1 .. II), j = 1 .. JJ)

H := 1/2*(1+tanh(w)):

plots:-display(plot3d(w, x = -1 .. 1, y = -1 .. 1, color = red), plot3d(H, x = -1 .. 1, y = -1 .. 1, color = blue)):

NULL

Warning,  computation interrupted

 

NULL

NULL

Download IntersectionIntegral.mw

I have black shaded the region as below figure that the integral must be computed inside it.

What do you suggest a fast and convenient approach to calculate a double integral inside the domain of intersection of two functions?

@tomleslie @tomleslie Let me explain my problem in a simpler way. I have a big 3-dimentional matrix (array) function of x[1], x[2], x[3], ..., x[n]. I am looking for the simplest command line way to transfer this matrix to Matlab. The small version of it seem as bellow:

restart;
A(1, 1, 1) := x[1]*x[2]*x[3];
A(2, 1, 1) := sin(x[1]*x[2]*x[3]);
A(1, 2, 1) := cos(x[1]*x[2]*x[3]);
A(2, 2, 1) := -x[2]*x[3]+x[1];
A(1, 1, 2) := x[1]*x[2]-x[3];
A(2, 1, 2) := x[1]-x[2]-x[3];
A(1, 2, 2) := x[1]-x[2]+x[3];
A(2, 2, 2) := -x[2]*x[3]+x[1];

Actually how is it possible to transfer a numeric and symbolic n-dimentional array from Maple to Matlab?

@tomleslie your approach looks nice, but I want to do it by command line not manually without the use of other softwares.

I used ExportMatrix, but faced error

ExportMatrix(matlabData, A, target = MATLAB, format = rectangular, mode = ascii);
Error, (in ExportMatrix) Matrix contains non-numeric entries:  cannot export to Matlab

@Carl Love Good comment, but not very usefull since this computation is only performed one time at the begining of the loop.

The main problem is that Maple allocates new memory for the same variables in each iteration without discarding the previous one. We need to find a way to discard this memory.

 

 

qs := proc (tn) local s; global t; Qn__2 := LinearAlgebra:-Transpose(Matrix(1, M+Nm)); Qn__3 := LinearAlgebra:-Transpose(Matrix(1, M+Nm)); F := diff((diff(sin(`&omega;__b`*t), t))*LinearAlgebra:-Transpose(convert(F0, Matrix)), t); s := 1; for t from t__0 by dt to tn do EQ := KL.V+eqMNL+a__0*ML.V+a__1*CL.V-ML.(a__0*q(1 .. M+Nm, s .. s)+a__2*dq(1 .. M+Nm, s .. s)+a__3*ddq(1 .. M+Nm, s .. s))-CL.(a__1*q(1 .. M+Nm, s .. s)+a__4*dq(1 .. M+Nm, s .. s)+a__5*ddq(1 .. M+Nm, s .. s))-F-Qn__2-Qn__3; G := Matrix(M+Nm, proc (i, j) options operator, arrow; (diff(EQ[i], V[j, 1]))[1] end proc); G1 := 1/G; if t = t__0 then X[1] := 1/(CL*a__1+ML*a__0+KL).(F+ML.(a__0*q(1 .. M+Nm, s .. s)+a__2*dq(1 .. M+Nm, s .. s)+a__3*ddq(1 .. M+Nm, s .. s))+CL.(a__1*q(1 .. M+Nm, s .. s)+a__4*dq(1 .. M+Nm, s .. s)+a__5*ddq(1 .. M+Nm, s .. s))) else X[1] := q(1 .. M+Nm, s .. s) end if; for k to 11 do X[k+1] := eval(V-G1.EQ, Equate(V, X[k])) end do; q(1 .. M+Nm, s+1 .. s+1) := X[k]; ddq(1 .. M+Nm, s+1 .. s+1) := a__0*(q(1 .. M+Nm, s+1 .. s+1)-q(1 .. M+Nm, s .. s))-a__2*dq(1 .. M+Nm, s .. s)-a__3*ddq(1 .. M+Nm, s .. s); dq(1 .. M+Nm, s+1 .. s+1) := dq(1 .. M+Nm, s .. s)+a__6*ddq(1 .. M+Nm, s .. s)+a__7*ddq(1 .. M+Nm, s+1 .. s+1); Wxy2 := add(add(add(h*W[i, j, 2, o]*LegendreP(i, `&zeta;__2`)*LegendreP(j, `&eta;__2`)*q[o, s+1], i = 0 .. II), j = 0 .. JJ), o = 1 .. M); Wxy3 := add(add(add(h*W[i, j, 3, o]*LegendreP(i, `&zeta;__2`)*LegendreP(j, `&eta;__2`)*q[o, s+1], i = 0 .. II), j = 0 .. JJ), o = 1 .. M); Plt := plots:-implicitplot(Wxy2-Wxy3, `&zeta;__2` = -1 .. 1, `&eta;__2` = -1 .. 1, color = red, thickness = 2, gridrefine = 3); j := 111111; Hvs2 := 1; for i while i <= j do Pnt[i] := op([1, i], Plt); if Wxy2-Wxy3 = 0 then j := 0 elif convert(Pnt[i], string)[1] = "C" then j := 0 else a__e[i] := (1/2)*abs(max(Column(Pnt[i], 1))-min(Column(Pnt[i], 1))); if a__e[i] = 0 then a__e[i] = 0.1e-5 end if; b__e[i] := (1/2)*abs(max(Column(Pnt[i], 2))-min(Column(Pnt[i], 2))); if b__e[i] = 0 then b__e[i] = 0.1e-5 end if; xc[i] := (1/2)*max(Column(Pnt[i], 1))+(1/2)*min(Column(Pnt[i], 1)); yc[i] := (1/2)*max(Column(Pnt[i], 2))+(1/2)*min(Column(Pnt[i], 2)); Hv2[i] := (`&zeta;__2`-xc[i])^2/a__e[i]^2+(`&eta;__2`-yc[i])^2/b__e[i]^2-1; Hvs2 := Hv2[i]*Hvs2 end if end do; Qnt := LinearAlgebra:-Transpose(convert([k__d*Grid:-Seq(Quadrature(Quadrature(Heaviside(-Hvs2)*abs(Wxy2-Wxy3)*add(add(W[i, j, 2, r]*LegendreP(i, `&zeta;__2`)*LegendreP(j, `&eta;__2`), i = 0 .. II), j = 0 .. JJ), `&zeta;__2` = -1 .. 1, method = romberg[5]), `&eta;__2` = -1 .. 1, method = romberg[5]), r = 1 .. M)/(a[2]*b[2]), seq(0, n = 1 .. Nm)], Matrix)); Qn__2 := Qnt; Qn__3 := -Qn__2; s := s+1; EQ := 'EQ'; G := 'G'; G1 := 'G1'; X[1] := 'X[1]'; Wxy2 := 'Wxy2'; Wxy3 := 'Wxy3'; Plt := 'Plt'; Hvs2 := 'Hvs2'; for ii to i do Pnt[ii] := 'Pnt[ii]'; a__e[ii] := 'a__e[ii]'; b__e[ii] := 'b__e[ii]'; xc[ii] := 'xc[ii]'; yc[ii] := 'yc[ii]'; Hv2[ii] := 'Hv2[ii]' end do; Hvs2 := 'Hvs'; Qnt := 'Qnt' end do; t := 't'; q end proc

NULL

NULL

with(CodeTools:-Profiling)

Profile(qs)

NULL

qs(dt)

PrintProfiles(qs)

qs
qs := proc(tn)
local s, Qn__2, Qn__3, F, EQ, G, G1, X, k, Wxy2, Wxy3, Plt, j, Hvs2, i, Pnt, a__e, b__e, xc, yc, Hv2, Qnt, ii;

global t;
     |Calls Seconds  Words|
PROC |    1   2.391 36561056|
   1 |    1   0.000    752| Qn__2 := LinearAlgebra:-Transpose(Matrix(1,M+Nm));
   2 |    1   0.000    752| Qn__3 := LinearAlgebra:-Transpose(Matrix(1,M+Nm));
   3 |    1   0.000   4009| F := MTM:-diff(MTM:-diff(MTM:-sin(omega__b*t),t)*LinearAlgebra:-Transpose(convert(F0,Matrix)),t);
   4 |    1   0.000      0| s := 1;
   5 |    1   0.000      4| for t from t__0 by dt to tn do
   6 |    1   0.016 180433|   EQ := Typesetting:-delayDotProduct(KL,V)+eqMNL+Typesetting:-delayDotProduct(a__0*ML,V)+Typesetting:-delayDotProduct(a__1*CL,V)-Typesetting:-delayDotProduct(ML,a__0*q(1 .. M+Nm,s .. s)+a__2*dq(1 .. M+Nm,s .. s)+a__3*ddq(1 .. M+Nm,s .. s))-Typesetting:-delayDotProduct(CL,a__1*q(1 .. M+Nm,s .. s)+a__4*dq(1 .. M+Nm,s .. s)+a__5*ddq(1 .. M+Nm,s .. s))-F-Qn__2-Qn__3;
   7 |    1   0.000   4376|   G := Matrix(M+Nm,(i, j) -> MTM:-diff(EQ[i],V[j,1])[1]);
   8 |    1   0.000   6755|   G1 := 1/G;
   9 |    1   0.000      3|   if t = t__0 then
  10 |    1   0.015 124077|     X[1] := Typesetting:-delayDotProduct(1/(CL*a__1+ML*a__0+KL),F+Typesetting:-delayDotProduct(ML,a__0*q(1 .. M+Nm,s .. s)+a__2*dq(1 .. M+Nm,s .. s)+a__3*ddq(1 .. M+Nm,s .. s))+Typesetting:-delayDotProduct(CL,a__1*q(1 .. M+Nm,s .. s)+a__4*dq(1 .. M+Nm,s .. s)+a__5*ddq(1 .. M+Nm,s .. s)))
                              else
  11 |    0   0.000      0|     X[1] := q(1 .. M+Nm,s .. s)
                              end if;
  12 |    1   0.000      0|   for k to 11 do
  13 |   11   0.016 376894|     X[k+1] := eval(V-Typesetting:-delayDotProduct(G1,EQ),Equate(V,X[k]))
                              end do;
  14 |    1   0.000     37|   q(1 .. M+Nm,s+1 .. s+1) := X[k];
  15 |    1   0.000   5925|   ddq(1 .. M+Nm,s+1 .. s+1) := a__0*(q(1 .. M+Nm,s+1 .. s+1)-q(1 .. M+Nm,s .. s))-a__2*dq(1 .. M+Nm,s .. s)-a__3*ddq(1 .. M+Nm,s .. s);
  16 |    1   0.000   5716|   dq(1 .. M+Nm,s+1 .. s+1) := dq(1 .. M+Nm,s .. s)+a__6*ddq(1 .. M+Nm,s .. s)+a__7*ddq(1 .. M+Nm,s+1 .. s+1);
  17 |    1   0.000 151003|   Wxy2 := add(add(add(h*W[i,j,2,o]*LegendreP(i,zeta__2)*LegendreP(j,eta__2)*q[o,s+1],i = 0 .. II),j = 0 .. JJ),o = 1 .. M);
  18 |    1   0.015 151003|   Wxy3 := add(add(add(h*W[i,j,3,o]*LegendreP(i,zeta__2)*LegendreP(j,eta__2)*q[o,s+1],i = 0 .. II),j = 0 .. JJ),o = 1 .. M);
  19 |    1   2.266 35071942|   Plt := plots:-implicitplot(Wxy2-Wxy3,zeta__2 = -1 .. 1,eta__2 = -1 .. 1,color = red,thickness = 2,gridrefine = 3);
  20 |    1   0.000      0|   j := 111111;
  21 |    1   0.000      0|   Hvs2 := 1;
  22 |    1   0.000     33|   for i while i <= j do
  23 |   10   0.000    433|     Pnt[i] := op([1, i],Plt);
  24 |   10   0.047  90747|     if Wxy2-Wxy3 = 0 then
  25 |    0   0.000      0|       j := 0
                                elif convert(Pnt[i],string)[1] = "C" then
  26 |    1   0.000      0|       j := 0
                                else
  27 |    9   0.000   7436|       a__e[i] := 1/2*MTM:-abs(max(LinearAlgebra:-Column(Pnt[i],1))-min(LinearAlgebra:-Column(Pnt[i],1)));
  28 |    9   0.000     81|       if a__e[i] = 0 then
  29 |    0   0.000      0|         a__e[i] = .1e-5
                                  end if;
  30 |    9   0.000   7113|       b__e[i] := 1/2*MTM:-abs(max(LinearAlgebra:-Column(Pnt[i],2))-min(LinearAlgebra:-Column(Pnt[i],2)));
  31 |    9   0.000     81|       if b__e[i] = 0 then
  32 |    0   0.000      0|         b__e[i] = .1e-5
                                  end if;
  33 |    9   0.000   7031|       xc[i] := 1/2*max(LinearAlgebra:-Column(Pnt[i],1))+1/2*min(LinearAlgebra:-Column(Pnt[i],1));
  34 |    9   0.000   6969|       yc[i] := 1/2*max(LinearAlgebra:-Column(Pnt[i],2))+1/2*min(LinearAlgebra:-Column(Pnt[i],2));
  35 |    9   0.000   1939|       Hv2[i] := (zeta__2-xc[i])^2/a__e[i]^2+(eta__2-yc[i])^2/b__e[i]^2-1;
  36 |    9   0.000    154|       Hvs2 := Hv2[i]*Hvs2
                                end if
                              end do;
  37 |    1   0.016 352769|   Qnt := LinearAlgebra:-Transpose(convert([k__d*Grid:-Seq(Student:-NumericalAnalysis:-Quadrature(Student:-NumericalAnalysis:-Quadrature(Heaviside(-Hvs2)*MTM:-abs(Wxy2-Wxy3)*add(add(W[i,j,2,r]*LegendreP(i,zeta__2)*LegendreP(j,eta__2),i = 0 .. II),j = 0 .. JJ),zeta__2 = -1 .. 1,method = romberg[5]),eta__2 = -1 .. 1,method = romberg[5]),r = 1 .. M)/a[2]/b[2], seq(0,n = 1 .. Nm)],Matrix));
  38 |    1   0.000      0|   Qn__2 := Qnt;
  39 |    1   0.000   1846|   Qn__3 := -Qn__2;
  40 |    1   0.000      0|   s := s+1;
  41 |    1   0.000      0|   EQ := 'EQ';
  42 |    1   0.000      0|   G := 'G';
  43 |    1   0.000      0|   G1 := 'G1';
  44 |    1   0.000      6|   X[1] := 'X[1]';
  45 |    1   0.000      0|   Wxy2 := 'Wxy2';
  46 |    1   0.000      0|   Wxy3 := 'Wxy3';
  47 |    1   0.000      0|   Plt := 'Plt';
  48 |    1   0.000      0|   Hvs2 := 'Hvs2';
  49 |    1   0.000      0|   for ii to i do
  50 |   11   0.000    117|     Pnt[ii] := 'Pnt[ii]';
  51 |   11   0.000    124|     a__e[ii] := 'a__e[ii]';
  52 |   11   0.000    124|     b__e[ii] := 'b__e[ii]';
  53 |   11   0.000    124|     xc[ii] := 'xc[ii]';
  54 |   11   0.000    124|     yc[ii] := 'yc[ii]';
  55 |   11   0.000    124|     Hv2[ii] := 'Hv2[ii]'
                              end do;
  56 |    1   0.000      0|   Hvs2 := 'Hvs';
  57 |    1   0.000      0|   Qnt := 'Qnt'
                            end do;
  58 |    1   0.000      0| t := 't';
  59 |    1   0.000      0| q
end proc

 

``

qs(2*dt)

PrintProfiles(qs)

qs
qs := proc(tn)
local s, Qn__2, Qn__3, F, EQ, G, G1, X, k, Wxy2, Wxy3, Plt, j, Hvs2, i, Pnt, a__e, b__e, xc, yc, Hv2, Qnt, ii;

global t;
     |Calls Seconds  Words|
PROC |    2   8.485 105246895|
   1 |    2   0.000   1504| Qn__2 := LinearAlgebra:-Transpose(Matrix(1,M+Nm));
   2 |    2   0.000   1504| Qn__3 := LinearAlgebra:-Transpose(Matrix(1,M+Nm));
   3 |    2   0.000   7448| F := MTM:-diff(MTM:-diff(MTM:-sin(omega__b*t),t)*LinearAlgebra:-Transpose(convert(F0,Matrix)),t);
   4 |    2   0.000      0| s := 1;
   5 |    2   0.000     12| for t from t__0 by dt to tn do
   6 |    3   0.016 543683|   EQ := Typesetting:-delayDotProduct(KL,V)+eqMNL+Typesetting:-delayDotProduct(a__0*ML,V)+Typesetting:-delayDotProduct(a__1*CL,V)-Typesetting:-delayDotProduct(ML,a__0*q(1 .. M+Nm,s .. s)+a__2*dq(1 .. M+Nm,s .. s)+a__3*ddq(1 .. M+Nm,s .. s))-Typesetting:-delayDotProduct(CL,a__1*q(1 .. M+Nm,s .. s)+a__4*dq(1 .. M+Nm,s .. s)+a__5*ddq(1 .. M+Nm,s .. s))-F-Qn__2-Qn__3;
   7 |    3   0.000  13136|   G := Matrix(M+Nm,(i, j) -> MTM:-diff(EQ[i],V[j,1])[1]);
   8 |    3   0.000  20270|   G1 := 1/G;
   9 |    3   0.000      9|   if t = t__0 then
  10 |    2   0.031 248154|     X[1] := Typesetting:-delayDotProduct(1/(CL*a__1+ML*a__0+KL),F+Typesetting:-delayDotProduct(ML,a__0*q(1 .. M+Nm,s .. s)+a__2*dq(1 .. M+Nm,s .. s)+a__3*ddq(1 .. M+Nm,s .. s))+Typesetting:-delayDotProduct(CL,a__1*q(1 .. M+Nm,s .. s)+a__4*dq(1 .. M+Nm,s .. s)+a__5*ddq(1 .. M+Nm,s .. s)))
                              else
  11 |    1   0.000     34|     X[1] := q(1 .. M+Nm,s .. s)
                              end if;
  12 |    3   0.000      0|   for k to 11 do
  13 |   33   0.078 1122449|     X[k+1] := eval(V-Typesetting:-delayDotProduct(G1,EQ),Equate(V,X[k]))
                              end do;
  14 |    3   0.000    111|   q(1 .. M+Nm,s+1 .. s+1) := X[k];
  15 |    3   0.000  17334|   ddq(1 .. M+Nm,s+1 .. s+1) := a__0*(q(1 .. M+Nm,s+1 .. s+1)-q(1 .. M+Nm,s .. s))-a__2*dq(1 .. M+Nm,s .. s)-a__3*ddq(1 .. M+Nm,s .. s);
  16 |    3   0.015  16801|   dq(1 .. M+Nm,s+1 .. s+1) := dq(1 .. M+Nm,s .. s)+a__6*ddq(1 .. M+Nm,s .. s)+a__7*ddq(1 .. M+Nm,s+1 .. s+1);
  17 |    3   0.016 453009|   Wxy2 := add(add(add(h*W[i,j,2,o]*LegendreP(i,zeta__2)*LegendreP(j,eta__2)*q[o,s+1],i = 0 .. II),j = 0 .. JJ),o = 1 .. M);
  18 |    3   0.047 453009|   Wxy3 := add(add(add(h*W[i,j,3,o]*LegendreP(i,zeta__2)*LegendreP(j,eta__2)*q[o,s+1],i = 0 .. II),j = 0 .. JJ),o = 1 .. M);
  19 |    3   8.031 101020050|   Plt := plots:-implicitplot(Wxy2-Wxy3,zeta__2 = -1 .. 1,eta__2 = -1 .. 1,color = red,thickness = 2,gridrefine = 3);
  20 |    3   0.000      0|   j := 111111;
  21 |    3   0.000      0|   Hvs2 := 1;
  22 |    3   0.000     99|   for i while i <= j do
  23 |   30   0.000    966|     Pnt[i] := op([1, i],Plt);
  24 |   30   0.141 271948|     if Wxy2-Wxy3 = 0 then
  25 |    0   0.000      0|       j := 0
                                elif convert(Pnt[i],string)[1] = "C" then
  26 |    3   0.000      0|       j := 0
                                else
  27 |   27   0.000  21307|       a__e[i] := 1/2*MTM:-abs(max(LinearAlgebra:-Column(Pnt[i],1))-min(LinearAlgebra:-Column(Pnt[i],1)));
  28 |   27   0.000    243|       if a__e[i] = 0 then
  29 |    0   0.000      0|         a__e[i] = .1e-5
                                  end if;
  30 |   27   0.000  21009|       b__e[i] := 1/2*MTM:-abs(max(LinearAlgebra:-Column(Pnt[i],2))-min(LinearAlgebra:-Column(Pnt[i],2)));
  31 |   27   0.000    243|       if b__e[i] = 0 then
  32 |    0   0.000      0|         b__e[i] = .1e-5
                                  end if;
  33 |   27   0.000  20597|       xc[i] := 1/2*max(LinearAlgebra:-Column(Pnt[i],1))+1/2*min(LinearAlgebra:-Column(Pnt[i],1));
  34 |   27   0.000  20577|       yc[i] := 1/2*max(LinearAlgebra:-Column(Pnt[i],2))+1/2*min(LinearAlgebra:-Column(Pnt[i],2));
  35 |   27   0.000   5480|       Hv2[i] := (zeta__2-xc[i])^2/a__e[i]^2+(eta__2-yc[i])^2/b__e[i]^2-1;
  36 |   27   0.000    462|       Hvs2 := Hv2[i]*Hvs2
                                end if
                              end do;
  37 |    3   0.110 957757|   Qnt := LinearAlgebra:-Transpose(convert([k__d*Grid:-Seq(Student:-NumericalAnalysis:-Quadrature(Student:-NumericalAnalysis:-Quadrature(Heaviside(-Hvs2)*MTM:-abs(Wxy2-Wxy3)*add(add(W[i,j,2,r]*LegendreP(i,zeta__2)*LegendreP(j,eta__2),i = 0 .. II),j = 0 .. JJ),zeta__2 = -1 .. 1,method = romberg[5]),eta__2 = -1 .. 1,method = romberg[5]),r = 1 .. M)/a[2]/b[2], seq(0,n = 1 .. Nm)],Matrix));
  38 |    3   0.000      0|   Qn__2 := Qnt;
  39 |    3   0.000   5538|   Qn__3 := -Qn__2;
  40 |    3   0.000      0|   s := s+1;
  41 |    3   0.000      0|   EQ := 'EQ';
  42 |    3   0.000      0|   G := 'G';
  43 |    3   0.000      0|   G1 := 'G1';
  44 |    3   0.000     18|   X[1] := 'X[1]';
  45 |    3   0.000      0|   Wxy2 := 'Wxy2';
  46 |    3   0.000      0|   Wxy3 := 'Wxy3';
  47 |    3   0.000      0|   Plt := 'Plt';
  48 |    3   0.000      0|   Hvs2 := 'Hvs2';
  49 |    3   0.000      0|   for ii to i do
  50 |   33   0.000    344|     Pnt[ii] := 'Pnt[ii]';
  51 |   33   0.000    358|     a__e[ii] := 'a__e[ii]';
  52 |   33   0.000    358|     b__e[ii] := 'b__e[ii]';
  53 |   33   0.000    358|     xc[ii] := 'xc[ii]';
  54 |   33   0.000    358|     yc[ii] := 'yc[ii]';
  55 |   33   0.000    358|     Hv2[ii] := 'Hv2[ii]'
                              end do;
  56 |    3   0.000      0|   Hvs2 := 'Hvs';
  57 |    3   0.000      0|   Qnt := 'Qnt'
                            end do;
  58 |    2   0.000      0| t := 't';
  59 |    2   0.000      0| q
end proc

 

``

``

Download Myqs.mw

@Joe Riel My main problem is more complicate, here I put my proc just to see how it is like. So it seems to me I only need to find a way for getting rid of new memmory allocation in loop.


 

``

qs := proc (tn) local s; global t; Qn__2 := LinearAlgebra:-Transpose(Matrix(1, M+Nm)); Qn__3 := LinearAlgebra:-Transpose(Matrix(1, M+Nm)); F := diff((diff(sin(`&omega;__b`*t), t))*LinearAlgebra:-Transpose(convert(F0, Matrix)), t); s := 1; for t from t__0 by dt to tn do EQ := KL.V+eqMNL+a__0*ML.V+a__1*CL.V-ML.(a__0*q(1 .. M+Nm, s .. s)+a__2*dq(1 .. M+Nm, s .. s)+a__3*ddq(1 .. M+Nm, s .. s))-CL.(a__1*q(1 .. M+Nm, s .. s)+a__4*dq(1 .. M+Nm, s .. s)+a__5*ddq(1 .. M+Nm, s .. s))-F-Qn__2-Qn__3; G := Matrix(M+Nm, proc (i, j) options operator, arrow; (diff(EQ[i], V[j, 1]))[1] end proc); G1 := 1/G; if t = t__0 then X[1] := 1/(CL*a__1+ML*a__0+KL).(F+ML.(a__0*q(1 .. M+Nm, s .. s)+a__2*dq(1 .. M+Nm, s .. s)+a__3*ddq(1 .. M+Nm, s .. s))+CL.(a__1*q(1 .. M+Nm, s .. s)+a__4*dq(1 .. M+Nm, s .. s)+a__5*ddq(1 .. M+Nm, s .. s))) else X[1] := q(1 .. M+Nm, s .. s) end if; for k to 11 do X[k+1] := eval(V-G1.EQ, Equate(V, X[k])) end do; q(1 .. M+Nm, s+1 .. s+1) := X[k]; ddq(1 .. M+Nm, s+1 .. s+1) := a__0*(q(1 .. M+Nm, s+1 .. s+1)-q(1 .. M+Nm, s .. s))-a__2*dq(1 .. M+Nm, s .. s)-a__3*ddq(1 .. M+Nm, s .. s); dq(1 .. M+Nm, s+1 .. s+1) := dq(1 .. M+Nm, s .. s)+a__6*ddq(1 .. M+Nm, s .. s)+a__7*ddq(1 .. M+Nm, s+1 .. s+1); Wxy2 := add(add(add(h*W[i, j, 2, o]*LegendreP(i, `&zeta;__2`)*LegendreP(j, `&eta;__2`)*q[o, s+1], i = 0 .. II), j = 0 .. JJ), o = 1 .. M); Wxy3 := add(add(add(h*W[i, j, 3, o]*LegendreP(i, `&zeta;__2`)*LegendreP(j, `&eta;__2`)*q[o, s+1], i = 0 .. II), j = 0 .. JJ), o = 1 .. M); Plt := plots:-implicitplot(Wxy2-Wxy3, `&zeta;__2` = -1 .. 1, `&eta;__2` = -1 .. 1, color = red, thickness = 2, gridrefine = 3); j := 111111; Hvs2 := 1; for i while i <= j do Pnt[i] := op([1, i], Plt); if Wxy2-Wxy3 = 0 then j := 0 elif convert(Pnt[i], string)[1] = "C" then j := 0 else a__e[i] := (1/2)*abs(max(Column(Pnt[i], 1))-min(Column(Pnt[i], 1))); if a__e[i] = 0 then a__e[i] = 0.1e-5 end if; b__e[i] := (1/2)*abs(max(Column(Pnt[i], 2))-min(Column(Pnt[i], 2))); if b__e[i] = 0 then b__e[i] = 0.1e-5 end if; xc[i] := (1/2)*max(Column(Pnt[i], 1))+(1/2)*min(Column(Pnt[i], 1)); yc[i] := (1/2)*max(Column(Pnt[i], 2))+(1/2)*min(Column(Pnt[i], 2)); Hv2[i] := (`&zeta;__2`-xc[i])^2/a__e[i]^2+(`&eta;__2`-yc[i])^2/b__e[i]^2-1; Hvs2 := Hv2[i]*Hvs2 end if end do; Qnt := LinearAlgebra:-Transpose(convert([k__d*Grid:-Seq(Quadrature(Quadrature(Heaviside(-Hvs2)*abs(Wxy2-Wxy3)*add(add(W[i, j, 2, r]*LegendreP(i, `&zeta;__2`)*LegendreP(j, `&eta;__2`), i = 0 .. II), j = 0 .. JJ), `&zeta;__2` = -1 .. 1, method = romberg[5]), `&eta;__2` = -1 .. 1, method = romberg[5]), r = 1 .. M)/(a[2]*b[2]), seq(0, n = 1 .. Nm)], Matrix)); Qn__2 := Qnt; Qn__3 := -Qn__2; s := s+1 end do; q end proc


 

Download MyProc.mw

I am wondering why Maple allocate new memory to the same variable in a loop without discarding the previous one.

@Joe Riel I created this analogy example since my main problem is something similar to it, but in a big scale. So your creative method does not work for my main problem. Any recommendation on getting rid of new memmory allocation?

@tomleslie as I emphasised in my question I am looking for numerical values of V[1], V[2], not symbolic. Your method only returns sympolic parametric values. Why the for loop in proc does not subs the values of t in equatins M:=A+A^(-1) and V[s]:=(M.B)???

@dharr I could not find helpful information in help page suitable for me, let me know if you have any sample Profiling or another suggestions.

@Carl Love I am wondering why the memory usage number at the the right side of the bottom is not compatible with the amount that task manager shows??!!!

Any change in new versions of Maple?

@acer Interesting idea.

@acer actually my function is composed of another factor as well, which by considering it your propose method does not work! What about proposing a proper numerical intergration method?

restart

kernelopts(version);

`Maple 2016.1, X86 64 WINDOWS, Apr 22 2016, Build ID 1133417`

(1)

A := Heaviside(zeta__2-1.000000000000)*Heaviside(eta__2+.9875792458758-3.881485663812*10^9*sqrt(9.765625000000*10^22-zeta__2^2))+Heaviside(zeta__2-.9999999999936)*Heaviside(eta__2+.9875792458758+3.881485663812*10^9*sqrt(9.765625000000*10^22-zeta__2^2))-Heaviside(zeta__2-.9999999999936)*Heaviside(eta__2+.9875792458758-3.881485663812*10^9*sqrt(9.765625000000*10^22-zeta__2^2))+Heaviside(zeta__2-.6466146460206)*Heaviside(eta__2-.5050000000000+8.127372424924*sqrt(269.5813999936-zeta__2^2))-Heaviside(zeta__2-.7684252323012)*Heaviside(eta__2-.5050000000000+8.127372424924*sqrt(269.5813999936-zeta__2^2))-Heaviside(zeta__2-.6466146460206)*Heaviside(eta__2-.5050000000000-8.127372424924*sqrt(269.5813999936-zeta__2^2))+Heaviside(zeta__2-.7684252323012)*Heaviside(eta__2-.5050000000000-8.127372424924*sqrt(269.5813999936-zeta__2^2))-Heaviside(zeta__2+1.000000000000)*Heaviside(eta__2+.7243706106403-1.382194833711*10^11*sqrt(1.562500000000*10^24-zeta__2^2))+Heaviside(zeta__2+1.000000000000)*Heaviside(eta__2+.9842650870048+1.085166413462*10^10*sqrt(4.756242568371*10^23-zeta__2^2))-Heaviside(zeta__2+1.000000000000)*Heaviside(eta__2+.9842650870048-1.085166413462*10^10*sqrt(4.756242568371*10^23-zeta__2^2))-Heaviside(zeta__2+.9999999999984)*Heaviside(eta__2-.9031048925918+8.123652892875*10^10*sqrt(1.562500000000*10^24-zeta__2^2))+Heaviside(zeta__2+.9999999999984)*Heaviside(eta__2-.9031048925918-8.123652892875*10^10*sqrt(1.562500000000*10^24-zeta__2^2))-Heaviside(zeta__2+.9999999999984)*Heaviside(eta__2+.7243706106403+1.382194833711*10^11*sqrt(1.562500000000*10^24-zeta__2^2))+Heaviside(zeta__2+.9999999999984)*Heaviside(eta__2+.7243706106403-1.382194833711*10^11*sqrt(1.562500000000*10^24-zeta__2^2))-Heaviside(zeta__2+.9999999999988)*Heaviside(eta__2-.4637698986762+2.327456686822*10^11*sqrt(2.777777777778*10^24-zeta__2^2))+Heaviside(zeta__2+.9999999999988)*Heaviside(eta__2-.4637698986762-2.327456686822*10^11*sqrt(2.777777777778*10^24-zeta__2^2))-Heaviside(zeta__2+.9999999999990)*Heaviside(eta__2+.1619291800251+3.018371923484*10^11*sqrt(4.000000000000*10^24-zeta__2^2))+Heaviside(zeta__2+.9999999999990)*Heaviside(eta__2+.1619291800251-3.018371923484*10^11*sqrt(4.000000000000*10^24-zeta__2^2))-Heaviside(zeta__2+.9999999999972)*Heaviside(eta__2+.9842650870048+1.085166413462*10^10*sqrt(4.756242568371*10^23-zeta__2^2))+Heaviside(zeta__2+.9999999999972)*Heaviside(eta__2+.9842650870048-1.085166413462*10^10*sqrt(4.756242568371*10^23-zeta__2^2))+Heaviside(zeta__2-.9999999999796)*Heaviside(eta__2-.8341191288491+1.298592148442*10^10*sqrt(9.706617486471*10^21-zeta__2^2))-Heaviside(zeta__2-.9999999999796)*Heaviside(eta__2-.8341191288491-1.298592148442*10^10*sqrt(9.706617486471*10^21-zeta__2^2))-Heaviside(zeta__2-1.000000000000)*Heaviside(eta__2-.8341191288491+1.298592148442*10^10*sqrt(9.706617486471*10^21-zeta__2^2))+Heaviside(zeta__2-1.000000000000)*Heaviside(eta__2-.8341191288491-1.298592148442*10^10*sqrt(9.706617486471*10^21-zeta__2^2))-Heaviside(zeta__2-1.000000000000)*Heaviside(eta__2+.9875792458758+3.881485663812*10^9*sqrt(9.765625000000*10^22-zeta__2^2))+Heaviside(zeta__2-.5527964251744)*Heaviside(eta__2+.5050000000000+10.98537767108*sqrt(492.5151416233-zeta__2^2))-Heaviside(zeta__2+1.000000000000)*Heaviside(eta__2+.1619291800251-3.018371923484*10^11*sqrt(4.000000000000*10^24-zeta__2^2))+Heaviside(zeta__2+1.000000000000)*Heaviside(eta__2+.7243706106403+1.382194833711*10^11*sqrt(1.562500000000*10^24-zeta__2^2))+Heaviside(zeta__2+1.000000000000)*Heaviside(eta__2-.4637698986762+2.327456686822*10^11*sqrt(2.777777777778*10^24-zeta__2^2))-Heaviside(zeta__2+1.000000000000)*Heaviside(eta__2-.4637698986762-2.327456686822*10^11*sqrt(2.777777777778*10^24-zeta__2^2))+Heaviside(zeta__2+1.000000000000)*Heaviside(eta__2+.1619291800251+3.018371923484*10^11*sqrt(4.000000000000*10^24-zeta__2^2))-Heaviside(zeta__2+1.000000000000)*Heaviside(eta__2-.9031048925918-8.123652892875*10^10*sqrt(1.562500000000*10^24-zeta__2^2))+Heaviside(zeta__2+1.000000000000)*Heaviside(eta__2-.9031048925918+8.123652892875*10^10*sqrt(1.562500000000*10^24-zeta__2^2))+Heaviside(zeta__2-.6429162216568)*Heaviside(eta__2+.5050000000000-10.98537767108*sqrt(492.5151416233-zeta__2^2))-Heaviside(zeta__2-.6429162216568)*Heaviside(eta__2+.5050000000000+10.98537767108*sqrt(492.5151416233-zeta__2^2))-Heaviside(zeta__2-.5527964251744)*Heaviside(eta__2+.5050000000000-10.98537767108*sqrt(492.5151416233-zeta__2^2)):

Digits := 22:

B := abs(1.54766124540895*10^(-23)*LegendreP(5, zeta__2)*eta__2+9.42848033882030*10^(-23)*LegendreP(4, zeta__2)*eta__2-8.63540720633871*10^(-25)*LegendreP(7, zeta__2)*eta__2+4.47592785073801*10^(-25)*LegendreP(6, zeta__2)*eta__2-3.27361791625699*10^(-24)*LegendreP(8, zeta__2)*eta__2+8.69201069226258*10^(-25)*LegendreP(9, zeta__2)*eta__2+8.44500836530084*10^(-24)*zeta__2*LegendreP(2, eta__2)-8.16447897662916*10^(-24)*LegendreP(8, zeta__2)-9.03325367966316*10^(-23)*LegendreP(2, zeta__2)*LegendreP(2, eta__2)-3.65665238354457*10^(-23)*LegendreP(3, zeta__2)*LegendreP(2, eta__2)+2.48616118186505*10^(-23)*LegendreP(5, zeta__2)*LegendreP(2, eta__2)+3.49274266801910*10^(-23)*LegendreP(4, zeta__2)*LegendreP(2, eta__2)+1.57562452147079*10^(-23)*LegendreP(7, zeta__2)*LegendreP(2, eta__2)+1.13836255333532*10^(-23)*LegendreP(6, zeta__2)*LegendreP(2, eta__2)-4.64840219234199*10^(-25)*LegendreP(9, zeta__2)*LegendreP(7, eta__2)-1.10434599924198*10^(-24)*zeta__2*LegendreP(8, eta__2)+1.08423343010855*10^(-23)*LegendreP(6, zeta__2)-1.80542010805784*10^(-23)*LegendreP(7, zeta__2)-2.49196685412280*10^(-24)*LegendreP(9, zeta__2)+2.43121828235083*10^(-22)*LegendreP(4, zeta__2)+8.00783752320630*10^(-25)*LegendreP(2, zeta__2)*LegendreP(3, eta__2)-5.24626201639296*10^(-25)*LegendreP(2, zeta__2)*LegendreP(8, eta__2)-4.46597278774917*10^(-25)*LegendreP(3, zeta__2)*LegendreP(8, eta__2)+4.40023693100841*10^(-22)*LegendreP(3, zeta__2)+3.96675680014788*10^(-25)*LegendreP(5, eta__2)-1.05763914360733*10^(-25)*LegendreP(8, eta__2)+1.01736979105438*10^(-24)*LegendreP(7, zeta__2)*LegendreP(3, eta__2)-2.11087035865237*10^(-24)*LegendreP(8, zeta__2)*LegendreP(3, eta__2)-6.55216211259737*10^(-24)*LegendreP(3, zeta__2)*LegendreP(4, eta__2)-2.59524263916176*10^(-24)*zeta__2*LegendreP(3, eta__2)-3.57670870238948*10^(-25)*LegendreP(2, zeta__2)*LegendreP(5, eta__2)+2.68358272136332*10^(-26)*LegendreP(7, eta__2)+9.25047899594019*10^(-25)*LegendreP(6, eta__2)-4.77479580653438*10^(-25)*zeta__2*LegendreP(5, eta__2)-2.71559028285720*10^(-26)*zeta__2*LegendreP(4, eta__2)+3.16345667898614*10^(-27)*LegendreP(9, eta__2)+1.86293670034105*10^(-22)*eta__2-1.39325020507283*10^(-24)*LegendreP(3, eta__2)-3.24357137541178*10^(-22)*zeta__2-4.54471957773919*10^(-24)*LegendreP(4, eta__2)+5.26383542096037*10^(-23)*LegendreP(2, eta__2)-7.87620119909178*10^(-22)*LegendreP(2, zeta__2)-9.51203876259311*10^(-23)*LegendreP(5, zeta__2)+7.99736443820554*10^(-25)*LegendreP(3, zeta__2)*LegendreP(6, eta__2)-1.71453452004739*10^(-25)*LegendreP(4, zeta__2)*LegendreP(6, eta__2)-2.18220802921933*10^(-24)*LegendreP(5, zeta__2)*LegendreP(6, eta__2)-3.77387134125069*10^(-24)*LegendreP(6, zeta__2)*LegendreP(6, eta__2)+1.00697041297227*10^(-24)*LegendreP(2, zeta__2)*LegendreP(6, eta__2)-2.92515100201635*10^(-24)*LegendreP(7, zeta__2)*LegendreP(6, eta__2)-3.68092156761947*10^(-24)*LegendreP(9, zeta__2)*LegendreP(4, eta__2)+1.74362117679247*10^(-24)*LegendreP(8, zeta__2)*LegendreP(4, eta__2)-3.32578513711331*10^(-25)*zeta__2*LegendreP(7, eta__2)+2.01330648090477*10^(-24)*LegendreP(8, zeta__2)*LegendreP(6, eta__2)+3.53233080368130*10^(-24)*LegendreP(9, zeta__2)*LegendreP(6, eta__2)-3.15725764921952*10^(-27)*LegendreP(4, zeta__2)*LegendreP(7, eta__2)+1.27534945335648*10^(-25)*LegendreP(2, zeta__2)*LegendreP(7, eta__2)-3.62425353476520*10^(-25)*LegendreP(3, zeta__2)*LegendreP(7, eta__2)+9.21474473171027*10^(-25)*LegendreP(7, zeta__2)*LegendreP(7, eta__2)+2.38369609942291*10^(-25)*LegendreP(5, zeta__2)*LegendreP(7, eta__2)-2.72256108440673*10^(-25)*LegendreP(6, zeta__2)*LegendreP(7, eta__2)-6.54138288634006*10^(-25)*LegendreP(7, zeta__2)*LegendreP(9, eta__2)+1.21042590275456*10^(-25)*LegendreP(8, zeta__2)*LegendreP(7, eta__2)+4.26720896569817*10^(-26)*LegendreP(6, zeta__2)*LegendreP(9, eta__2)-5.12858575529453*10^(-26)*LegendreP(8, zeta__2)*LegendreP(9, eta__2)+4.76504374040493*10^(-25)*LegendreP(9, zeta__2)*LegendreP(9, eta__2)+3.37897478439262*10^(-23)*zeta__2*eta__2-2.77752448289485*10^(-22)*LegendreP(2, zeta__2)*eta__2-3.08054377508128*10^(-24)*LegendreP(6, zeta__2)*LegendreP(4, eta__2)+1.15853047423759*10^(-23)*LegendreP(7, zeta__2)*LegendreP(4, eta__2)-4.92720206449747*10^(-23)*LegendreP(3, zeta__2)*eta__2-8.61686962535103*10^(-24)*LegendreP(8, zeta__2)*LegendreP(2, eta__2)-1.32506515971502*10^(-24)*LegendreP(5, zeta__2)*LegendreP(4, eta__2)+6.85899626186726*10^(-25)*LegendreP(5, zeta__2)*LegendreP(3, eta__2)+2.28203090969654*10^(-24)*LegendreP(3, zeta__2)*LegendreP(3, eta__2)-2.68907301809976*10^(-25)*LegendreP(5, zeta__2)*LegendreP(9, eta__2)+3.63485480297361*10^(-26)*LegendreP(4, zeta__2)*LegendreP(9, eta__2)+1.97348199257558*10^(-25)*LegendreP(3, zeta__2)*LegendreP(9, eta__2)-3.08982363240662*10^(-26)*LegendreP(2, zeta__2)*LegendreP(9, eta__2)+2.49193016043161*10^(-25)*zeta__2*LegendreP(9, eta__2)-4.18880639432511*10^(-24)*LegendreP(9, zeta__2)*LegendreP(8, eta__2)-2.96092073149334*10^(-25)*LegendreP(8, zeta__2)*LegendreP(8, eta__2)-1.39005768761275*10^(-24)*LegendreP(9, zeta__2)*LegendreP(3, eta__2)+7.75291783959110*10^(-25)*zeta__2*LegendreP(6, eta__2)+8.64278679290519*10^(-25)*LegendreP(6, zeta__2)*LegendreP(8, eta__2)+3.54095069410550*10^(-24)*LegendreP(7, zeta__2)*LegendreP(8, eta__2)+2.19879897667752*10^(-24)*LegendreP(5, zeta__2)*LegendreP(8, eta__2)+6.22035083218788*10^(-26)*LegendreP(4, zeta__2)*LegendreP(8, eta__2)+2.72240385770454*10^(-25)*LegendreP(8, zeta__2)*LegendreP(5, eta__2)-1.93048147157364*10^(-24)*LegendreP(4, zeta__2)*LegendreP(3, eta__2)-1.24963415620315*10^(-23)*LegendreP(9, zeta__2)*LegendreP(2, eta__2)+4.63381828314138*10^(-24)*LegendreP(6, zeta__2)*LegendreP(3, eta__2)+1.02854711982595*10^(-26)*LegendreP(4, zeta__2)*LegendreP(5, eta__2)-1.08656057300493*10^(-24)*LegendreP(3, zeta__2)*LegendreP(5, eta__2)+1.01562215037002*10^(-23)*LegendreP(2, zeta__2)*LegendreP(4, eta__2)+2.22669013695569*10^(-24)*LegendreP(7, zeta__2)*LegendreP(5, eta__2)-3.21530667993533*10^(-25)*LegendreP(6, zeta__2)*LegendreP(5, eta__2)+1.07214426103800*10^(-24)*LegendreP(5, zeta__2)*LegendreP(5, eta__2)+5.41820436348695*10^(-22)-4.27457932806028*10^(-24)*LegendreP(4, zeta__2)*LegendreP(4, eta__2)-1.73479424553310*10^(-24)*LegendreP(9, zeta__2)*LegendreP(5, eta__2)):

plot3d(A*B, `&zeta;__2` = -1 .. 1, `&eta;__2` = -1 .. 1, color = green)

with(Student[NumericalAnalysis]):

Quadrature(Quadrature(A*B, zeta__2 = -1 .. 1, method = romberg[3]), eta__2 = -1 .. 1, method = romberg[3])

Float(undefined)

(2)

forget(`evalf/int`);

Warning,  computation interrupted

 

forget(`evalf/int`); forget(evalf); CodeTools:-Usage(evalf(Int(A*B, [`&zeta;__2` = -1 .. 1, `&eta;__2` = -1 .. 1], method = _CubaCuhre, epsilon = 0.1e-3))); evalf[4](%)

 

Download romberg_ac(1).mw

@dharr The appeared example is a simple case study. If fact, there is no guarantee for linearity of equations.

1 2 3 4 5 6 7 Last Page 1 of 9