|
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]);
|
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;
|
The Eigenvalues and Eigenvectors computed with MATLAB® are found:
|
>
|
(P,W):=eig(K,M,eigenvectors=true):
|
The Eigenvalues are:
The Eigenvectors are:
The equivalent computations in the Maple environment:
|
>
|
(W2,P2) := Eigenvectors(K,M):
|
The Eigenvalues are:
the Eigenvectors are (in no particular order):
|
|
|
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;
|
|
>
|
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] );
|
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");
|
Or, we can solve the system entirely within Maple:
|
>
|
U2 := MatrixInverse(A) . B;
|
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(kernelopts(homedir));
|
|
>
|
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))) ];"):
|
|
>
|
writeline(file, "end"):
|
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]);
|
|
>
|
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)]:
|
|
>
|
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):
|
|
>
|
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");
|
|
>
|
evalM("x=sin(2*pi*50*t) + sin(2*pi*120*t)");
|
|
>
|
evalM("y=x+2*randn(size(t))");
|
|
>
|
evalM("Pyy=Y.*conj(Y)/512");
|
|
>
|
evalM("f=1000*(0:511)/512");
|
|
>
|
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
|
