|
Maple-MATLAB Link Application
Initiate the MATLAB link. Use the with command to access the functions in some useful packages by their short names:
|
>
|
restart:
with(LinearAlgebra):
with(plots):
with(Matlab):
|
For more information on the Maple-MATLAB link, see Matlab.
|
|
Structural Analysis: A First Approximation
|
|
This formulation is used to compute the lowest natural frequencies and modes of a highly idealized 22-story building.
The model equations may be formulated with the following matrix assignments.
The mass matrix M is a matrix with m on the diagonal:
|
>
|
m := 5000.:
M := DiagonalMatrix([seq(m,i=1..22)], outputoptions=[datatype=float,shape=symmetric,attributes=[positive_definite],storage=rectangular]);
|
|
| (1.1) |
To examine any of these Matrices, the Structured Data Browser can be used, by right-clicking the output Matrix and selecting Browse.
The Stiffness matrix K is tridiagonal with 2k on the center diagonal, and -k on the adjacent diagonals:
|
>
|
k := 1.25e8:
K := BandMatrix([-k,2*k,-k],1,22,22, outputoptions=[datatype=float,shape=symmetric,storage=rectangular]):
K[22,22] := -k:
K;
|
|
| (1.2) |
The Eigenvalues and Eigenvectors computed with Matlab are found:
|
>
|
(P,W):=eig(K,M,eigenvectors=true):
|
The Eigenvalues are:
|
| (1.3) |
The Eigenvectors are:
|
| (1.4) |
The equivalent computations in the Maple environment:
|
>
|
(W2,P2) := Eigenvectors(K,M):
|
The Eigenvalues are:
|
| (1.5) |
the Eigenvectors are (in no particular order):
|
| (1.6) |
|
|
|
Heat Transfer: Finite Difference Solution
|
|
A difference equation method is used to find the static temperature distribution in a flat rectangular plate, given its boundary is held at a fixed temperature profile.
We model the plate as a 3 x 7 grid of nodes, where the nodes may be thought of as being interconnected with a square mesh of heat conductors.
With our 21 plate-internal nodes and 20 boundary conditions ( = 7+3+7+3 ), the finite-difference 2-dimensional Laplace equation gives us an inhomogeneous linear system in 21 unknowns (the internal nodal temperatures). Representing the internal nodal temperatures U[i,j] by the column vector [ U[1,1], ..., U[1,7], U[2,1], ..., U[2,7], U[3,1], ..., U[3,7] ],we have the system AU=B, where the matrix A and column vector B encode the interconnectivity of the nodes and the profile of the boundary temperature:
|
>
|
A := BandMatrix([1,0,0,0,0,0,1,-4,1,0,0,0,0,0,1],7,21,21,outputoptions=[datatype=float[8]]):
A[7,8]:=0:A[8,7]:=0:
A[14,15]:=0:A[15,14]:=0:
A;
|
|
| (2.1) |
|
>
|
T := 100.:
B := Vector[column]( 21, [-T,-T,-T,-T,-T,-T,-T,0,0,0,0,0,0,0,0,0,0,0,0,0,0], datatype=float[8] );
|
|
| (2.2) |
We have fixed the temperature along the U[1,1], ..., U[1,7] edge at 100 units, and the other three edges at 0 temperature units.
Now, we solve the system in the MATLAB environment:
|
>
|
evalM("U=inv(A) * B");
U1 := getvar("U");
|
|
| (2.3) |
Or, we can solve the system entirely within Maple:
|
>
|
U2 := MatrixInverse(A) . B;
|
|
| (2.4) |
Now, let us take the Maple solution vector and re-arrange its values in a 3 x 7 matrix that corresponds to the actual rectangular plate:
|
>
|
P := <<U2[1..7]>|<U2[8..14]>|<U2[15..21]>>:
|
|
>
|
plots[matrixplot](P,heights=histogram,axes=boxed,labels=["","","T"],title="temperature distribution");
|
|
|
|
Simulation of a Linear Oscillator
|
|
A simple spring-mass-dashpot is modeled as a second-order linear oscillator. The simulation equations are coded as MATLAB function files and are called from the Maple environment by using the ode45 command.
The MATLAB function stored in the file mass_eqn.m can be created within Maple as follows, or it can be created by using any text editor. Here, we create the file in the current directory.
|
>
|
currentdir("/usr/local");
|
|
>
|
file := open("mass_eqn.m", WRITE):
|
|
>
|
writeline(file, "function xdot=mass_eqn(t,x)"):
writeline(file, " global M C K"):
writeline(file, " xdot = [ x(2); (1/M*(-K*x(1)-C*x(2))) ];"):
|
The oscillator parameters of Mass M, Damping C, and Stiffness K are:
Pass the variables to MATLAB. We must enforce the above variables as global in the MATLAB environment so that the function mass_eqn can use them.
|
>
|
setvar("M", M, 'globalvar'):
setvar("C", C, 'globalvar'):
setvar("K", K, 'globalvar'):
|
Allow MATLAB to solve the ODE.
|
>
|
(t,x) := ode45("mass_eqn", 0..50, [2,0]);
|
|
| (3.2) |
|
>
|
plots[pointplot]([seq([t[i],x[i,1]],i=1..Dimensions(t))], symbol=DIAMOND, symbolsize=8);
|
|
|
|
Data Analysis: Fast Fourier Transform (FFT)
|
|
A set of experimental numerical data can be analyzed with MATLAB as a numerical engine.
|
>
|
t := [seq(i*.001, i=0..600)]:
|
|
| (4.1) |
|
>
|
x := map(proc(T) evalhf(sin(2*Pi*50*T) + sin(2*Pi*120*T)) end, t):
|
|
>
|
data := map(proc(X) X+2*stats[random,normald](1) end, x):
|
|
| (4.2) |
|
>
|
Pyy := [seq(abs(Y[i])^2/ 512, i=1..512)]:
|
|
>
|
freq := [seq(1000*i/512., i=0..511)]:
|
|
>
|
plots[pointplot]([seq([t[i],data[i]],i=1..nops(t))],symbol=DIAMOND,symbolsize=7);
|
|
>
|
plots[pointplot]([seq([freq[i],Pyy[i]],i=1..512)],style=line);
|
Alternatively, the same calculations can be done by using MATLAB syntax almost entirely.
|
>
|
evalM("t=0:0.001:0.6");
|
|
| (4.3) |
|
| (4.4) |
|
>
|
evalM("x=sin(2*pi*50*t) + sin(2*pi*120*t)");
|
|
>
|
evalM("y=x+2*randn(size(t))");
|
|
| (4.5) |
|
>
|
evalM("Pyy=Y.*conj(Y)/512");
|
|
| (4.6) |
|
>
|
evalM("f=1000*(0:511)/512");
|
|
| (4.7) |
|
| (4.8) |
|
| (4.9) |
|
>
|
plots[pointplot]([seq([t[i],data[i]],i=1..tmax)],symbol=DIAMOND,symbolsize=7);
|
|
>
|
plots[pointplot]([seq([freq[i],Pyy[i]],i=1..512)],style=line);
|
|
Return to Index for Example Worksheets
|
