# Using APDLMath to solve Eigenproblems#

Use APDLMath to solve eigenproblems.

This example uses a verification manual input file, but you can use your own sparse or dense matrices and solve those.

import time

import matplotlib.pylab as plt
import numpy as np

from ansys.mapdl.core import launch_mapdl
from ansys.mapdl.core.examples import vmfiles

# Start MAPDL as a service and create an APDLMath object
mapdl = launch_mapdl(loglevel="ERROR")
mm = mapdl.math


First we get the STIFF and MASS matrices from the full file after running the input file from Verification Manual 153

out = mapdl.input(vmfiles["vm153"])

k = mm.stiff(fname="PRSMEMB.full")
m = mm.mass(fname="PRSMEMB.full")


Display size of the M and K matrices

print(m.shape)
print(k.shape)


Out:

(126, 126)
(126, 126)


Allocate an array to store the eigenshapes. where nev is the number of eigenvalues requested

nev = 10
a = mm.mat(k.nrow, nev)
a


Out:

Dense APDLMath Matrix (126, 10)


Perform the the modal analysis.

The algorithm is automatically chosen with respect to the matrices properties (e.g. scalar, storage, symmetry…)

print("Calling MAPDL to solve the eigenproblem...")

t1 = time.time()
ev = mm.eigs(nev, k, m, phi=a)
print(f"Elapsed time to solve this problem: {time.time() - t1}")


Out:

Calling MAPDL to solve the eigenproblem...
/opt/hostedtoolcache/Python/3.8.12/x64/lib/python3.8/site-packages/ansys/mapdl/core/math.py:1453: UserWarning: Call to sym cannot evaluate if this matrix is symmetric with this version of MAPDL.
warn(
Elapsed time to solve this problem: 0.030767202377319336


This is the vector of eigenfrequencies.

print(ev)


Out:

FQAMOZ :
Size : 10
3.381e+02   3.381e+02   6.266e+02   6.266e+02   9.283e+02      <       5
9.283e+02   1.250e+03   1.250e+03   1.424e+03   1.424e+03      <       10


## Verify the accuracy of eigenresults#

Check the residual error for the first eigenresult $$R_1=||(K-\lambda_1.M).\phi_1||_2$$

First, we compute $$\lambda_1 = \omega_1^2 = (2.\pi.f_1)^2$$

# Eigenfrequency (Hz)
i = 0
f = ev[0]
omega = 2 * np.pi * f
lam = omega * omega


Then we get the 1st Eigenshape $$\phi_1$$, and compute $$K.\phi_1$$ and $$M.\phi_1$$

# shape
phi = a[0]

# APDL Command: *MULT,K,,Phi,,KPhi
kphi = k.dot(phi)

# APDL Command: *MULT,M,,Phi,,MPhi
mphi = m.dot(phi)


Next, compute the $$||K.\phi_1||_2$$ quantity and normalize the residual value.

# APDL Command: *MULT,K,,Phi,,KPhi
kphi = k.dot(phi)

# APDL Command: *NRM,KPhi,NRM2,KPhiNrm
kphinrm = kphi.norm()


Then we add these two vectors, using the $$\lambda_1$$ scalar factor and finally compute the normalized residual value $$\frac{R_1}{||K.\phi_1||_2}$$

# APDL Command: *AXPY,-lambda,,MPhi,1,,KPhi
mphi *= lam
kphi -= mphi

# Compute the residual
res = kphi.norm() / kphinrm
print(res)


Out:

7.11275926845182e-11


This residual can be computed for all eigenmodes

def get_res(i):
"""Compute the residual for a given eigenmode"""
# Eigenfrequency (Hz)
f = ev[i]

# omega = 2.pi.Frequency
omega = 2 * np.pi * f

# lambda = omega^2
lam = omega * omega

# i-th eigenshape
phi = a[i]

# K.Phi
kphi = k.dot(phi)

# M.Phi
mphi = m.dot(phi)

# Normalization scalar value
kphinrm = kphi.norm()

# (K-\lambda.M).Phi
mphi *= lam
kphi -= mphi

# return the residual
return kphi.norm() / kphinrm

mapdl_acc = np.zeros(nev)

for i in range(nev):
f = ev[i]
mapdl_acc[i] = get_res(i)
print(f"[{i}] : Freq = {f}\t - Residual = {mapdl_acc[i]}")


Out:

[0] : Freq = 338.06666355063646  - Residual = 7.11275926845182e-11
[1] : Freq = 338.06666355063675  - Residual = 6.746171206219663e-11
[2] : Freq = 626.645098092703    - Residual = 1.7213067209485443e-11
[3] : Freq = 626.6450980927034   - Residual = 3.4196464654380337e-11
[4] : Freq = 928.2598500574524   - Residual = 6.447818437869047e-12
[5] : Freq = 928.2598500574531   - Residual = 1.2195745273859991e-11
[6] : Freq = 1249.8421074363494  - Residual = 1.9669862802830984e-11
[7] : Freq = 1249.8421074363505  - Residual = 1.8811190184372674e-11
[8] : Freq = 1423.993890941667   - Residual = 3.132913660044634e-10
[9] : Freq = 1423.9938909416703  - Residual = 1.2587507950775475e-09


Plot Accuracy of Eigenresults

fig = plt.figure(figsize=(12, 10))
ax = plt.axes()
x = np.linspace(1, nev, nev)
plt.title("APDL Math Residual Error (%)")
plt.yscale("log")
plt.ylim([10e-13, 10e-7])
plt.xlabel("Frequency #")
plt.ylabel("Errors (%)")
ax.bar(x, mapdl_acc, label="MAPDL Results")


Out:

<BarContainer object of 10 artists>


stop mapdl

mapdl.exit()


Total running time of the script: ( 0 minutes 1.090 seconds)

Gallery generated by Sphinx-Gallery