Maple Questions and Posts

These are Posts and Questions associated with the product, Maple

As an addition to the post.
Non-orientable surface in the sequence of orientable surfaces. In the picture we see the equations corresponding to the current surface plot.
Just entertainment.
surfaces.mw

Hello,

I have a question about fitting a function to experiemental data. The experimental data is Di = X(zi) +I*Y(zi) for a zi = [z1 z1+dz ...z1+dz*N], where zi are real values. I want to fit a complex function f(z,a,b,c) to the data zi and Di , where a,b,c are complex-valued parameters. Does someone know how to fit a complex function to a complex data? The Maple function that I found work only with real functions and data.

Hello,

I want to define an orthonormal tetrad basis of my choice in a spacetime having a metric given in some system of coordinates. My problem is that Maple automatically proposes an orthonormal metric but this is not the one that suits my requirements. So, I would like to specify the tetrad basis manually. As an example, I am trying to reproduce the calculations in sections 6 and 7 of the article https://arxiv.org/abs/gr-qc/0510083 . Here, the metric $g$ is given by the line element $ds^2 = - (c(t,r)^2 - v(t,r)^2) dt^2 + 2 v(t,r) dr dt + dr^2 + r^2 (d\theta^2 + sin(\theta)^2 d\phi^2)$ in $(t, r, \theta, \phi)$ coordinates. My chosen signature is (- + + +). Let, us adopt the convention used by Maple and denote spacetime indices by Greek alphabets and tetrad indices by lowercase Latin letters. Now, I would like to define a tetrad $e_a = (V, S, \Theta, \Phi)$ (as in section 7 of the article referred to above) where:

V^\mu = \frac{1}{c\sqrt{1-\beta(t,r)^2}}[1, - (v + c \beta), 0, 0] \\

S^\mu = \frac{1}{c\sqrt{1-\beta^2}}[-\beta, c + v \beta, 0, 0] \\

\Theta^\mu = [0,0,1,0]

\Phi^\mu = [0,0,0,1].

Here, $|\beta(t,r)| < 1$. I do not know how I may specify this in my worksheet. This may come of use somewhere later. Now, with this choice of the tetrad, we know that $g(e_a, e_b) = \eta_{ab}$ with $\eta$ being the Minkowski metric in spherical coordinates. After defining this tetrad basis, I finally want to calculate Einstein tensor, components of energy-momentum tensr etc. I have problem with constructing this orthonormal tetrad basis myself. It would be great if you could help me with this.

 

An additional curiosity: when we work with multiple tetrad bases, is it possible to denote the the tetrad indices by hatted tetrad labels themselves, as in $\eta_{\hat V, \hat \Theta}$?

 

Thank you.
 

restart

with(Physics)

[`*`, `.`, Annihilation, AntiCommutator, Antisymmetrize, Assume, Bra, Bracket, Cactus, Check, Christoffel, Coefficients, Commutator, CompactDisplay, Coordinates, Creation, D_, Dagger, Decompose, Define, Dgamma, Einstein, EnergyMomentum, Expand, ExteriorDerivative, Factor, FeynmanDiagrams, Fundiff, Geodesics, GrassmannParity, Gtaylor, Intc, Inverse, Ket, KillingVectors, KroneckerDelta, LeviCivita, Library, LieBracket, LieDerivative, Normal, Parameters, PerformOnAnticommutativeSystem, Projector, Psigma, Redefine, Ricci, Riemann, Setup, Simplify, SpaceTimeVector, StandardModel, SubstituteTensor, SubstituteTensorIndices, SumOverRepeatedIndices, Symmetrize, TensorArray, Tetrads, ThreePlusOne, ToFieldComponents, ToSuperfields, Trace, TransformCoordinates, Vectors, Weyl, `^`, dAlembertian, d_, diff, g_, gamma_]

(1)

Setup(signature = `-+++`, coordinates = (X = [t, r, theta, phi]))

`* Partial match of  'coordinates' against keyword 'coordinatesystems'`

 

`Default differentiation variables for d_, D_ and dAlembertian are: `*{X = (t, r, theta, phi)}

 

`Systems of spacetime Coordinates are: `*{X = (t, r, theta, phi)}

 

[coordinatesystems = {X}, signature = `- + + +`]

(2)

Setup(g_=-(c(t,r)^2 - v(t,r)^2)*dt^2 + 2*v(t,r)*dt*dr + dr^2 + r^2*dtheta^2 + r^2*sin(theta)^2*dphi^2)

[metric = {(1, 1) = -c(t, r)^2+v(t, r)^2, (1, 2) = v(t, r), (2, 2) = 1, (3, 3) = r^2, (4, 4) = r^2*sin(theta)^2}]

(3)

PDETools:-declare(c(t, r), v(t, r))

` c`(t, r)*`will now be displayed as`*c

 

` v`(t, r)*`will now be displayed as`*v

(4)

with(Tetrads)

`Setting lowercaselatin_ah letters to represent tetrad indices `

 

0, "%1 is not a command in the %2 package", Tetrads, Physics

 

0, "%1 is not a command in the %2 package", Tetrads, Physics

 

[IsTetrad, NullTetrad, OrthonormalTetrad, PetrovType, SegreType, TransformTetrad, e_, eta_, gamma_, l_, lambda_, m_, mb_, n_]

(5)

e_[]

Physics:-Tetrads:-e_[a, mu] = Matrix(%id = 18446744078438692614)

(6)

``


 

Download dynBH.mw

 


 

 

 

#include<iostream>
#include<vector>
#include<cmath>
#define NODE 8

using namespace std;
int graph[NODE][NODE] = {
   {0,1,1,0,0,0,0,0},
   {1,0,1,1,1,0,0,0},
   {1,1,0,1,0,1,0,0},
   {0,1,1,0,0,0,0,0},
   {0,1,0,0,0,1,1,1},
   {0,0,1,0,1,0,1,1},
   {0,0,0,0,1,1,0,0},
   {0,0,0,0,1,1,0,0}
};
int tempGraph[NODE][NODE];
int findStartVert() {
   for(int i = 0; i<NODE; i++) {
      int deg = 0;
      for(int j = 0; j<NODE; j++) {
         if(tempGraph[i][j])
            deg++; //increase degree, when connected edge found
      }
      if(deg % 2 != 0) //when degree of vertices are odd
      return i; //i is node with odd degree
   }
   return 0; //when all vertices have even degree, start from 0
}
int dfs(int prev, int start, bool visited[]){
   int count = 1;
   visited[start] = true;
   for(int u = 0; u<NODE; u++){
      if(prev != u){
         if(!visited[u]){
            if(tempGraph[start][u]){
               count += dfs(start, u, visited);
            }
         }
      }
   }
   return count;
}
bool isBridge(int u, int v) {
   int deg = 0;
   for(int i = 0; i<NODE; i++)
      if(tempGraph[v][i])
   deg++;
   if(deg>1) {
      return false; //the edge is not forming bridge
   }
   return true; //edge forming a bridge
}
int edgeCount() {
   int count = 0;
   for(int i = 0; i<NODE; i++)
      for(int j = i; j<NODE; j++)
         if(tempGraph[i][j])
   count++;
   return count;
}
void fleuryAlgorithm(int start) {
   static int edge = edgeCount();
   static int v_count = NODE;
   for(int v = 0; v<NODE; v++) {
      if(tempGraph[start][v]) {
         bool visited[NODE] = {false};
         if(isBridge(start, v)){
            v_count--;
         }
         int cnt = dfs(start, v, visited);
         if(abs(v_count-cnt) <= 2){
            cout << start << "--" << v << " ";
            if(isBridge(v, start)){
               v_count--;
            }
            tempGraph[start][v] = tempGraph[v][start] = 0; //remove edge from graph
            edge--;
            fleuryAlgorithm(v);
         }
      }
   }
}
int main() {
   for(int i = 0; i<NODE; i++) //copy main graph to tempGraph
   for(int j = 0; j<NODE; j++)
      tempGraph[i][j] = graph[i][j];
   cout << "Euler Path Or Circuit: ";
   fleuryAlgorithm(findStartVert());
}

Kind help 

 

 

Hi, I generated latex formate of an equation by using a command of maple but when I paste it into MathType, could not get the required equation, can anyone help me

${\frac {1}{51200\, \left( {x}^{2}+2 \right) ^{6}} \left( -187110\,

 \left( {x}^{2}+2 \right) ^{6}\sqrt {2} \left( {Q}^{3}+ \left( {\frac

{18\,k}{11}}-{\frac{18}{11}} \right) {Q}^{2}+ \left( {\frac {320\,{k}^

{2}}{297}}-{\frac {40\,k}{27}}+{\frac{320}{297}} \right) Q+{\frac {80

\,{k}^{3}}{297}}-{\frac {80\,{k}^{2}}{189}}+{\frac {80\,k}{189}}+{

\frac {640\,\lambda}{2079}}-{\frac{80}{297}} \right) \arctan \left( 1/

2\,x\sqrt {2} \right) -93555\, \left( {x}^{2}+2 \right) ^{6}\pi\,

 \left( {Q}^{3}+ \left( {\frac {18\,k}{11}}-{\frac{18}{11}} \right) {Q

}^{2}+ \left( {\frac {320\,{k}^{2}}{297}}-{\frac {40\,k}{27}}+{\frac{

320}{297}} \right) Q+{\frac {80\,{k}^{3}}{297}}-{\frac {80\,{k}^{2}}{

189}}+{\frac {80\,k}{189}}+{\frac {640\,\lambda}{2079}}-{\frac{80}{297

}} \right) \sqrt {2}-374220\, \left(  \left( {Q}^{3}+ \left( {\frac {

18\,k}{11}}-{\frac{18}{11}} \right) {Q}^{2}+ \left( {\frac {320\,{k}^{

2}}{297}}-{\frac {40\,k}{27}}+{\frac{320}{297}} \right) Q+{\frac {80\,

{k}^{3}}{297}}-{\frac {80\,{k}^{2}}{189}}+{\frac {80\,k}{189}}+{\frac

{640\,\lambda}{2079}}-{\frac{80}{297}} \right) {x}^{10}+ \left( {

\frac {34\,{Q}^{3}}{3}}+ \left( {\frac {204\,k}{11}}-{\frac{204}{11}}

 \right) {Q}^{2}+ \left( {\frac {10880\,{k}^{2}}{891}}-{\frac {1360\,k

}{81}}+{\frac{10880}{891}} \right) Q+{\frac {2720\,{k}^{3}}{891}}-{

\frac {2720\,{k}^{2}}{567}}+{\frac {2720\,k}{567}}+{\frac {21760\,

\lambda}{6237}}-{\frac{2720}{891}} \right) {x}^{8}+ \left( {\frac {264

\,{Q}^{3}}{5}}+ \left( {\frac {432\,k}{5}}-{\frac{432}{5}} \right) {Q}

^{2}+ \left( {\frac {512\,{k}^{2}}{9}}-{\frac {704\,k}{9}}+{\frac{512}

{9}} \right) Q+{\frac {128\,{k}^{3}}{9}}-{\frac {1408\,{k}^{2}}{63}}+{

\frac {1408\,k}{63}}+{\frac {97280\,\lambda}{6237}}-{\frac{128}{9}}

 \right) {x}^{6}+ \left( {\frac {4496\,{Q}^{3}}{35}}+ \left( {\frac {

80928\,k}{385}}-{\frac{80928}{385}} \right) {Q}^{2}+ \left( {\frac {

287744\,{k}^{2}}{2079}}-{\frac {35968\,k}{189}}+{\frac{287744}{2079}}

 \right) Q+{\frac {3328\,{k}^{3}}{99}}-{\frac {3328\,{k}^{2}}{63}}+{

\frac {3328\,k}{63}}+{\frac {10240\,\lambda}{297}}-{\frac{3328}{99}}

 \right) {x}^{4}+ \left( {\frac {10672\,{Q}^{3}}{63}}+ \left( {\frac {

21344\,k}{77}}-{\frac{21344}{77}} \right) {Q}^{2}+ \left( {\frac {

1094656\,{k}^{2}}{6237}}-{\frac {136832\,k}{567}}+{\frac{1094656}{6237

}} \right) Q+{\frac {35584\,{k}^{3}}{891}}-{\frac {35584\,{k}^{2}}{567

}}+{\frac {35584\,k}{567}}+{\frac {235520\,\lambda}{6237}}-{\frac{

35584}{891}} \right) {x}^{2}+{\frac {25376\,{Q}^{3}}{231}}+ \left( {

\frac {12352\,k}{77}}-{\frac{12352}{77}} \right) {Q}^{2}+ \left( -{

\frac {7936\,k}{63}}+{\frac {63488\,{k}^{2}}{693}}+{\frac{63488}{693}}

 \right) Q-{\frac{512}{27}}+{\frac {512\,{k}^{3}}{27}}-{\frac {5632\,{

k}^{2}}{189}}+{\frac {102400\,\lambda}{6237}}+{\frac {5632\,k}{189}}

 \right) x \right) }$

 

 

 

Hi,

I am trying to build a procedure to check to see if an nxn matrix is elementary. I am working with the standard definition of an elementary matrix, a square matrix that is formed by performing exactly one row operation on the corresponding identity matrix. It could also be interpreted as a matrix that is exactly one row operation from becoming the identity matrix. Does anyone have any idea how this might be done?

Thanks!

Hello Everyone,

I am currently in the process of building a user interface where a user can input a ply layup for a composite part which needs to be calculated. I want to forbid the use of non-numerical value in the matrix which represents my ply layup, before it gets sent to be calculated, but I am wondering how to check if my matrix has non-numerical values in it. 

At best, it would not only check if the matrix has non-numerical entries, but also give me the indices where exactly that value is present. I have yet been unsuccessful trying to find a fitting command. Do you perhaps know something that might work.
For clairification, allowable entries are in a range of -90 to 90, so it should not look for the "-" in the matrix. 

Best Regards,
Lennart

Hello Everyone,

So I am currently plotting some data within an embedded plot component using the following code:

SetProperty("DeterminantenRaumPlot1", value, dataplot(DetAnalysisNx2, DetAnalysisDeterm2));

Which is working just fine, DetAnalysisNx2 and DetAnalysisDeterm2 are single rows vectors with numbers in it.

However I need the plot to be in logarithmic mode for the vertical axis, but I do not know how to achieve this. I can set it manually, but when I save and reopen the document and re-do the calculation it turns back to linear vertical axis.

What kind of code segment am I missing? I have tried adding this like 'mode = log' but it is not working with dataplot. Maybe I need a different kind of plot? But then my Vectors will not work anymore, dataplot makes it so easy to use and the setting is available manually.

 

Any help would be fantastic!
Thanks and best regards,
Lennart

Hello all, 

Would you allow me to ask this question?

Is there a way to make the expression 'expression1', given below, to expression 'expression3'?


 

restart;

expression1 := -cos(theta)/2 + sin(theta)*sqrt(3)/2;

-(1/2)*cos(theta)+(1/2)*sin(theta)*3^(1/2)

(1)

expression2 := expand(cos(theta-2/3*Pi));

-(1/2)*cos(theta)+(1/2)*sin(theta)*3^(1/2)

(2)

is(expression1 = expression2);

true

(3)

expression3 := cos(theta-2/3*Pi);

-cos(theta+(1/3)*Pi)

(4)

 


Thank you, 

Download Q20201007.mw

I have four equations,

y = 1/x, y=3, y=1, x=0

restart;

P1 := plots:-implicitplot(y=1/x, x=-1..2, y=0..4, color=black, scaling=constrained):
P2 := plots:-implicitplot(y=3,   x=-1..2, y=0..4, color=black):
P3 := plots:-implicitplot(y=1,   x=-1..2, y=0..4, color=black):

plots:-display([P1,P2,P3]);

How do I fill the area between the first three equations and the y axis (x=0)?

Is it possible to do this with plottools:-transform?

How to obtain the multiple solution and graph given in the paper. 

Stefan Blowing and Slip Effects on Unsteady Nanofluid Transport Past a Shrinking Sheet: Multiple Solutions

https://doi.org/10.1002/htj.21470

 

Can anyone help to get solutions.

We recently published a paper on multiscale-multidomain simulation of battery models.

https://iopscience.iop.org/article/10.1149/1945-7111/abb37b

Some challenges are listed at 

https://twitter.com/UT_MAPLE/status/1311356957941522438

Maple and symbolic math can play a critical role in solving many challenging problems. For example, consider a seemingly-trivial problem

uxx+uyy = 0

x = 0 and x = 1, ux = 0 for all y

y = 1, u = 0, for all x

y = 0, u =1, 0<=x<=0.5

y = 0, uy = 0, 0.5<x<=1

There is a singularity at (0.5,0) and most numerical methods will have trouble there. In these equations, uxx means the second derivative of u with respect to x.

Maple can help solve this problem with conformal mapping to achieve arbitrary precision. As of today, machine precision is not possible with any numerical method even with millions of Degrees of Freedom. The Maple code is given below. A FEM code is given below as well. Models like this can benefit from Maple adding parallel sparse direct and iterative solvers.


 

 

restart;

Digits:=15;

Digits := 15

(1)

 The domain is tranformed from Z = X+IY to w. The points tranformed are

Zdomain:=[[0,0],[1,0],[1,1],[0,1]];

Zdomain := [[0, 0], [1, 0], [1, 1], [0, 1]]

(2)

wdomain:=[0,1,1+a,a+2];

wdomain := [0, 1, 1+a, a+2]

(3)

eq1:=diff(Z(w),w)=-I*K1/sqrt(w)/sqrt(w-1)/sqrt(w-a-1)/sqrt(w-a-2);

eq1 := diff(Z(w), w) = -I*K1/(w^(1/2)*(w-1)^(1/2)*(w-a-1)^(1/2)*(w-a-2)^(1/2))

(4)

 a is not known apriori. The value of a should be found to make sure [1,1] in the Z coordinate is transformed to [1+a] in the w coordinate.

a:=sqrt(2)-1;

a := 2^(1/2)-1

(5)

 Value of K1 is found using the transformation of [1,0] to 1 in the w coordinate

eq11:=1=int(rhs(eq1),w=0..1);

eq11 := 1 = -K1*2^(1/2)*EllipticK(2^(1/2)/2)

(6)

K1:=solve(eq11,K1);

K1 := -(1/2)*(2^(1/2)/EllipticK(2^(1/2)/2))

(7)

 The height Y in the Z coordinate is found by integrating from 1 to 1+a in the w coordinate.

simplify(int(rhs(eq1),w=1..1+a));

(2*I)*2^(1/2)*EllipticK((2^(1/2)-1)/(2^(1/2)+1))/(EllipticK(2^(1/2)/2)*(2^(1/2)+1))

(8)

evalf(%);

1.00000000000000*I

(9)

 The choice of a = sqrt(2)-1 gives the height of 1 for Z coordinate.

eval(simplify(int(rhs(eq1),w=0..1+a)));

(2^(1/2)*EllipticK(2^(1/2)/2)+EllipticK(2^(1/2)/2)+(2*I)*2^(1/2)*EllipticK((2^(1/2)-1)/(2^(1/2)+1)))/(EllipticK(2^(1/2)/2)*(2^(1/2)+1))

(10)

evalf(%);

1.00000000000000+1.00000000000001*I

(11)

 integration from 0 to 1+a in the w coordinate gives 1,1 in the Z coordinate.

simplify(int(rhs(eq1),w=0..1/sqrt(2)));

EllipticF(2^(1/2)/(2+2^(1/2))^(1/2), 2^(1/2)/2)/EllipticK(2^(1/2)/2)

(12)

evalf(%);

.500000000000001

(13)

 Integrating w from 0 to wmid =1/sqrt(2) gives the point 0.5,0 in the Z coordinate

wmid:=1/sqrt(2);

wmid := 2^(1/2)/2

(14)

 Next w domain is transformed to Z2 domain Z2 = X2+IY2

Z2domain:=[[0,0],[1,0],[1,H],[0,H]];

Z2domain := [[0, 0], [1, 0], [1, H], [0, H]]

(15)

wdomain2:=[0,1/sqrt(2),1+a,a+2];

wdomain2 := [0, 2^(1/2)/2, 2^(1/2), 2^(1/2)+1]

(16)

eq2:=diff(Z2(w),w)=-I*K2/sqrt(w)/sqrt(w-wmid)/sqrt(w-1-a)/sqrt(w-a-2);

eq2 := diff(Z2(w), w) = -(2*I)*K2/(w^(1/2)*(4*w-2*2^(1/2))^(1/2)*(w-2^(1/2))^(1/2)*(w-2^(1/2)-1)^(1/2))

(17)

 K2 is found based on the transformation of wmid to [1,0] in Z2 coordinate.

eq21:=1=int(rhs(eq2),w=0..wmid);

eq21 := 1 = -2*K2*EllipticK(1/(((2+2^(1/2))^(1/2))))/(2^(1/2)+1)^(1/2)

(18)

K2:=solve(eq21,K2);

K2 := -(1/2)*((2^(1/2)+1)^(1/2)/EllipticK(1/(((2+2^(1/2))^(1/2)))))

(19)

 The corner 0,1 in the Z coordinate is mapped by integrating eq2 from 0 to 1 in the w coordinate

 

 

int(rhs(eq2),w=0..1);

1+EllipticF((2-2^(1/2))^(1/2), ((2^(1/2)+1)/(2+2^(1/2)))^(1/2))*I/EllipticK(1/(((2+2^(1/2))^(1/2))))

(20)

corner:=evalf(%);

corner := 1.+.559419351518322*I

(21)

 This is the point 1,0 in the original coordinate.

 The height in the Z2 coorinate is found by integrating eq2 from wmid to 1+a.

ytot:=int(rhs(eq2),w=wmid..1+a);

ytot := EllipticK(((2^(1/2)+1)/(2+2^(1/2)))^(1/2))*I/EllipticK(1/(((2+2^(1/2))^(1/2))))

(22)

evalf(%);

1.22004159128347*I

(23)

 The magnitude in the Y direction is given by the coefficient of I, the imaginary number

Ytot:=coeff(ytot,I);

Ytot := EllipticK(((2^(1/2)+1)/(2+2^(1/2)))^(1/2))/EllipticK(1/(((2+2^(1/2))^(1/2))))

(24)

 The analytial solution in the Z2 corordinate is a line in Y2 to satisfy simple zero flux conditions at X2 = 0 and X2 = 1.

phianal:=1+b*Y2;

phianal := 1+b*Y2

(25)

 The value of phi is zero at Y2 = Ytot (originally the cathode domain in the Z domain)

bc:=subs(Y2=Ytot,phianal)=0;

bc := 1+b*EllipticK(((2^(1/2)+1)/(2+2^(1/2)))^(1/2))/EllipticK(1/(((2+2^(1/2))^(1/2)))) = 0

(26)

b:=solve(bc,b);

b := -EllipticK(1/(((2+2^(1/2))^(1/2))))/EllipticK(((2^(1/2)+1)/(2+2^(1/2)))^(1/2))

(27)

 The analytical solution is given by

phianal;

1-EllipticK(1/(((2+2^(1/2))^(1/2))))*Y2/EllipticK(((2^(1/2)+1)/(2+2^(1/2)))^(1/2))

(28)

evalf(phianal);

1.-.819644188480505*Y2

(29)

 The potential at the corner is given by substituting the imaginary value of corner for Y2 in phinanal)

phicorner:=subs(Y2=Im(corner),phianal);

phicorner := 1-.559419351518322*EllipticK(1/(((2+2^(1/2))^(1/2))))/EllipticK(((2^(1/2)+1)/(2+2^(1/2)))^(1/2))

(30)

evalf(phicorner);

.541475179604475

(31)

local flux/current density calculation, written in terms of w is

curr:=b*rhs(eq2)/rhs(eq1);

curr := -(2^(1/2)+1)^(1/2)*2^(1/2)*EllipticK(2^(1/2)/2)*(w-1)^(1/2)/(EllipticK(((2^(1/2)+1)/(2+2^(1/2)))^(1/2))*(4*w-2*2^(1/2))^(1/2))

(32)

average flux/current density calculation for the anode

currave:=int((curr),w=0..wmid)/wmid;

currave := -(1/2)*((2^(1/2)+1)^(1/2)*(2^(1/2)-1)*EllipticK(2^(1/2)/2)*(2^(3/4)+2^(1/4)+arctanh(2^(3/4)/2))*2^(1/2)/EllipticK(2^(3/4)/2))

(33)

Digits:=25:

 The average current density at Y =0, local current density at X = 0,Y=0 and potential at X=1,Y=0 (Corner) can be used to study convergence of FEM and other numerical methods

evalf(currave),evalf(subs(w=0,curr)),evalf(phicorner);

-1.656507648777793388522396, -1.161311530233258689567781, .5414751796044734741869534

(34)

 


 

Download Conformalmapping.mws


 

This FEM code is for solving Laplace's equation with primary current distribution considered in Model 1.
This code is based on FEM weak-form. Biquadratic Lagrange shape functions (9nodes in an element) are used.

restart;

with(LinearAlgebra):

Lx:=1: #length in X

Ly:=1: #length in Y

nx:=10: #number of elements in X (even numbers only)

ny:=10: #number of elements in Y, to be kept same as nx in this version

hx:=Lx/nx: #element size x

hy:=Ly/ny: #element size y

Procedure to perform numerical integration on shape functions to obtain local matrices (can be replaced with analytical integration for this particular problem)
  -Shape functions are also used as weight functions in applying weak formulation. Numerical integration is done using Simpson's rule.
  -Local cartesian coordinates x,y are converted to natural coordinates zeta and eta. This transformation is not required for this simple geomerty but useful in general. zeta and eta are obtained by scaling x and y with hx/2 and hy/2, respectively, in this code.

 

localmatrices:=proc(a1,a2,a3,b1,b2,b3,q1,q2)
global Kx,Ky,Nx,Ny,zeta,eta,c;
local A,dAdzeta,dAdeta,y,x,J,terms,i,j,k,l,dx,dy,fx,fy,fxy,fyy,dzeta,deta,J1,J2;

A:=[(1-zeta)*zeta*(1-eta)*eta/4,-(1-(zeta)^2)*(1-eta)*eta/2,-(1+zeta)*zeta*(1-eta)*eta/4,(1-(eta)^2)*(1+zeta)*zeta/2,(1+zeta)*zeta*(1+eta)*eta/4,(1-(zeta)^2)*(1+eta)*eta/2,-(1-zeta)*zeta*(1+eta)*eta/4,-(1-(eta)^2)*(1-zeta)*(zeta)/2,(1-(zeta)^2)*(1-(eta)^2)]; #bi quadratic langrange shape functions

dAdzeta:=diff(A,zeta);

dAdeta:=diff(A,eta);

y:=[-a1,-a2,-a3,0,a3,a2,a1,0,0];x:=[-b1,0,b1,b2,b3,0,-b3,-b2,0];

J:=simplify(add(dAdzeta[i]*x[i],i=1..9)*add(dAdeta[i]*y[i],i=1..9)-add(dAdzeta[i]*y[i],i=1..9)*add(dAdeta[i]*x[i],i=1..9));

Nx:=[seq((simplify(add(dAdeta[i]*y[i],i=1..9))*dAdzeta[j]-simplify(add(dAdzeta[i]*y[i],i=1..9))*dAdeta[j])/simplify(J),j=1..9)];

Ny:=[seq((-dAdzeta[j]*simplify(add(dAdeta[i]*x[i],i=1..9))+dAdeta[j]*simplify(add(dAdzeta[i]*x[i],i=1..9)))/simplify(J),j=1..9)];

Kx:=Matrix(nops(A),nops(A),datatype=float[8]):

Ky:=Matrix(nops(A),nops(A),datatype=float[8]):
c:=Vector(nops(A),datatype=float[8]):

terms:=20:#number of terms for numerical integration
dzeta:=2/terms:
deta:=2/terms:

for i from 1 to nops(Nx) do #loop to obtain local matrices      

for j from 1 to nops(Ny) do
Kx[i,j]:=0;
Ky[i,j]:=0;

for k from 0 to terms do #outer loop double integration, integration in zeta

if k = 0 then fx[k]:= subs(zeta=-1,Nx[i]*Nx[j]*J); fy[k]:= subs(zeta=-1,Ny[i]*Ny[j]*J);  
elif k = terms then
fx[k]:= subs(zeta=-1+(k*dzeta),Nx[i]*Nx[j]*J);
fy[k]:= subs(zeta=-1+(k*dzeta),Ny[i]*Ny[j]*J);  
elif irem(k,2) = 0 then
fx[k]:= 2*subs(zeta=-1+(k*dzeta),Nx[i]*Nx[j]*J);
fy[k]:=     2*subs(zeta=-1+(k*dzeta),Ny[i]*Ny[j]*J);
else fx[k]:= 4*subs(zeta=-1+(k*dzeta),Nx[i]*Nx[j]*J);
fy[k]:=     4*subs(zeta=-1+(k*dzeta),Ny[i]*Ny[j]*J);  end if;

for l from 0 to terms do #inner loop double integration, integration in eta

if l = 0 then fxy[l]:= subs(eta=-1,fx[k]); fyy[l]:= subs(eta=-1,fy[k]);
elif l = terms then fxy[l]:= subs(eta=-1+(l*deta),fx[k]); fyy[l]:= subs(eta=-1+(l*deta),fy[k]);
elif irem(l,2) = 0 then fxy[l]:= 2*subs(eta=-1+(l*deta),fx[k]); fyy[l]:= 2*subs(eta=-1+(l*deta),fy[k]);
else fxy[l]:=4*subs(eta=-1+(l*deta),fx[k]); fyy[l]:=4*subs(eta=-1+(l*deta),fy[k]); end if;
Kx[i,j]:=Kx[i,j]+fxy[l];
Ky[i,j]:=Ky[i,j]+fyy[l];

end do;
    
end do;
Kx[i,j]:=Kx[i,j]*dzeta*deta/9;
Ky[i,j]:=Ky[i,j]*dzeta*deta/9;
end do;
end do:

end proc:

n:=nx*ny; #total number of elements

n := 100

(1)

Nx1:=nx*2+1: #number of nodes in x in one row

N:=Nx1*(2*ny+1); # total number of nodes/equations

N := 441

(2)

K:=Matrix(N,N,storage=sparse): # global K matrix

C:=Vector(N,storage=sparse): # global c matrix

L2G:=Matrix(n,9):  #mapping matrix - each row has node numbers for each element

l:=1:k:=1:
localmatrices(hy/2,hy/2,hy/2,hx/2,hx/2,hx/2,0,0): kx:=copy(Kx):ky:=copy(Ky):c0:=copy(c):

for i from 1 to n do #modifying,adding and assembling matrices to get global matrix
 
if i<=nx/2 then  
a1:=copy(kx); a2:=copy(ky); a3:=0; a1[1..3,1..9]:=IdentityMatrix(3,9); a2[1..3,1..9]:=Matrix(3,9,shape=zero); a4:=a1+a2; c:=copy(c0); c[1..3]:=1.0;

elif i=nx/2+1 then  a1:=copy(kx); a2:=copy(ky); a3:=0; a1[1,1..9]:=IdentityMatrix(1,9); a2[1,1..9]:=Matrix(1,9,shape=zero); a4:=a1+a2; c:=copy(c0); c[1]:=1.0;

elif i>nx*(ny-1) then a1:=copy(kx); a2:=copy(ky); a3:=0; a1[5..7,5..7]:=IdentityMatrix(3,3);
a1[5..7,1..4]:=ZeroMatrix(3,4);
a1[5..7,8..9]:=ZeroMatrix(3,2); a2[5..7,1..9]:=Matrix(3,9,shape=zero); a4:=a1+a2; c:=copy(c0); c[5..7]:=0;

else a1:=kx; a2:=ky; a3:=0;a4:=a1+a2; c:=c0;  end if;

L2G[i,1..9]:=Matrix([l,l+1,l+2,l+2+Nx1,l+2+Nx1*2,l+1+Nx1*2,l+Nx1*2,l+Nx1,l+1+Nx1]):

k:=k+1:
 
if k>nx then k:=1; l:=l+Nx1+3;

else l:=l+2; end if:

indx2:=L2G[i,1..9]:
indx2:=convert(indx2,list):

C[indx2]:=C[indx2]+c[1..9];
c[1..9];

for i1 from 1 to 9 do
indx1:=L2G[i,i1]:
K[indx1,indx2]:=K[indx1,indx2]+a4[i1,1..9]:
end do:

end do:

phi:=LinearSolve(K,C,method=SparseLU): #linear set of equations solved using Sparse LU solver

phi_at(1,0):=phi[Nx1];
 

phi_at(1, 0) := .546587799122513

(3)

dNdy:=copy(Ny):

dNdy_bottom_left:=subs(eta=-1,zeta=-1,dNdy):

current_at(0,0):=add(dNdy_bottom_left[i]*phi[L2G[1,i]],i=1..nops(Ny));

current_at(0, 0) := -1.15773815354626

(4)

if irem(nx/2,2)=0
then current_at(0,0.25):=add(dNdy_bottom_left[i]*phi[L2G[nx/4+1,i]],i=1..nops(Ny));
else
dNdy_bottom_center:=subs(eta=-1,zeta=0,dNdy):
current_at(0,0.25):=add(dNdy_bottom_center[i]*phi[L2G[(nx/2+1)/2,i]],i=1..nops(Ny));
end if;

dNdy_bottom_center := [0, -30, 0, 0, 0, -10, 0, 0, 40]

current_at(0, .25) := -1.26989097821724

(5)

 

 


 

Download FEM_2D.mws

Hi,

It might be a really silly question, but I am wondering is it possible to simplify expression like this

a^2+b^2+2*a*b+c^2. Just by looking it we know that we can write it in the form of (a+b)^2+c^2. This is the basic exmaple I come up with. I have very lengthy expressions in maple, which can be factorize like this, but factor command will not work as it will try to factorize entire expression. So I am wondering is it doable in maple or I have to do it manually by collecting terms and check whether they can be factorize or not.

Thanks in Advance.

With Regards

sunit

I know this has to do with name scoping issue. But I do not see how to fix it now.

Calling J,Q:=LinearAlgebra:-JordanForm(A,output=['J','Q']); works in global, but not inside a proc.

What is the correct way to use this inside a proc?


 

restart;
A:=Matrix([[1,2],[3,4]]);
J,Q:=LinearAlgebra:-JordanForm(A,output=['J','Q']);

Matrix(2, 2, {(1, 1) = 1, (1, 2) = 2, (2, 1) = 3, (2, 2) = 4})

Matrix(%id = 18446745500542083790), Matrix(%id = 18446745500542076678)

restart;
foo:=proc(A::Matrix)
local J,Q;
J,Q:=LinearAlgebra:-JordanForm(A,output=['J','Q']);
end proc;
A:=Matrix([[1,2],[3,4]]);
foo(A);

foo := proc (A::Matrix) local J, Q; J, Q := LinearAlgebra:-JordanForm(A, output = ['J', 'Q']) end proc

Matrix(%id = 18446745500542658910)

Error, (in foo) invalid input: LinearAlgebra:-JordanForm expects value for keyword parameter output to be of type {list(identical(J, Q)), identical(J, Q)}, but received [J, Q]

 


 

Download jordan_issue.mw

Hi,

I would like to instal Maple 9.5 in my laptop, once I already have some few programs for his version. I don't know how much cost Maple 9.5, and how do download it. 

I live in Brazil.

Thanks in advance.

Carrijo.

jose.carrijo@gmail.com

2 3 4 5 6 7 8 Last Page 4 of 1737