第1章：MATLAB编程语言概述
- 二分法编程示例
% 二分法求根示例
function root = bisection_method(a, b, tol)
% 检查区间两端点的函数值是否具有相反的符号
if sign(func(a)) == sign(func(b))
error('区间两端点的函数值必须具有相反的符号。');
end
% 迭代直到达到指定的容差
while (b - a) / 2 > tol
c = (a + b) / 2;
% 检查中点是否已经是根
if func(c) == 0
root = c;
return;
elseif sign(func(c)) == sign(func(a))
a = c;
else
b = c;
end
end
root = (a + b) / 2;
end
% 定义要求根的函数
function y = func(x)
y = x^3 - 2*x - 5;
end
- 牛顿-拉夫森法编程示例
% 牛顿-拉夫森法求根示例
function root = newton_raphson_method(x0, tol)
% 迭代直到达到指定的容差
while true
fx = func(x0);
dfx = deriv_func(x0);
% 计算下一个近似值
x1 = x0 - fx / dfx;
% 检查是否满足收敛条件
if abs(x1 - x0) < tol
root = x1;
return;
end
x0 = x1;
end
end
% 定义要求根的函数
function y = func(x)
y = x^3 - 2*x - 5;
end
% 定义要求根函数的导数
function y = deriv_func(x)
y = 3*x^2 - 2;
end
- 高斯-牛顿算法编程示例（非线性回归）
% 高斯-牛顿算法用于非线性回归示例
% Kian Aghani & Salar Farahmand-Tabar (2023)
clear; clc
% 输入数据
y = [0.05 0.127 0.094 0.2122 0.2729 0.2665 0.3317];
x = [0.038 0.194 0.426 0.626 1.253 2.5 3.74];
m = 2;
n = length(x);
a0 = ones(2, 1);
maxiter = 500;
jacobian = zeros(n, m);
errorMat = zeros(n, 1);
toler = 0.0001;
for iter = 1:maxiter
for k = 1:n
dRda = -x(k) / (x(k) + a0(2));
dTdb = x(k) * a0(1) / (x(k) + a0(2))^2;
jacobian(k, 1) = dRda;
jacobian(k, 2) = dTdb;
Ytest = x(k) * a0(1) / (x(k) + a0(2));
errorMat(k) = y(k) - Ytest;
end
normError = norm(errorMat);
delta = - (jacobian' * jacobian) \ jacobian' * errorMat;
a0 = a0 + delta;
plot(iter, normError, 'k--o');
hold on;
if norm(delta) < toler
h = ['收敛在：', num2str(iter), '次迭代后'];
disp(h);
break;
end
end
xlabel('迭代');
ylabel('误差');
disp(a0);
fun = @(a0, x) a0(1) * x./ (a0(2) + x);
figure;
plot(x, y, 'bo', x, fun(a0, x), 'k');
legend('数据', '回归');
xlabel('x');
ylabel('y');
第2章：框架结构的矩阵分析
- 平面桁架相关函数
function [S] = planetruss(E, A, Wd, element, elementcon, xcord, ycord)
S = zeros(Wd, Wd);
for i = 1:element
dof = elementcon(i, :);
x = xcord(dof(2)) - xcord(dof(1));
y = ycord(dof(2)) - ycord(dof(1));
L = sqrt(x^2 + y^2);
l = x / L;
m = y / L;
R = [l m 0 0; 0 0 l m];
k = (E * A / L) * [1 -1; -1 1];
K = R' * k * R;
index = [dof(1) * 2 - 1 dof(1) * 2 dof(2) * 2 - 1 dof(2) * 2];
S(index, index) = S(index, index) + K;
end
end
- `drawplanetruss`函数
function drawplanetruss(element, elementcon, nodecord, newcord, Snewcord)
a = max(Snewcord(:, 1));
b = min(nodecord(:, 1));
x0 = b(1) - 2;
x1 = a(1) * 1.2;
y0 = 0;
y1 = a(2) * 1.2;
for n = 1:element
dof = elementcon(n, :);
xx = [nodecord(dof(1), 1), nodecord(dof(2), 1)];
yy = [nodecord(dof(1), 2), nodecord(dof(2), 2)];
xxn = [Snewcord(dof(1), 1), Snewcord(dof(2), 1)];
yyn = [Snewcord(dof(1), 2), Snewcord(dof(2), 2)];
plot(xx, yy, '-.k', xxn, yyn, 'k');
grid on;
title('Scaled Deformation');
xlabel('Horizontal Deformation');
ylabel('Vertical deformation');
axis([x0, x1, y0, y1]);
hold on;
end
end
- `trussforce`函数
function [p] = trussforce(E, A, element, elementcon, xcord, ycord, U)
p = zeros(element, 1);
for j = 1:element
dof = elementcon(j, :);
x = xcord(dof(2)) - xcord(dof(1));
y = ycord(dof(2)) - ycord(dof(1));
L = sqrt(x^2 + y^2);
l = x / L;
m = y / L;
index = [dof(1) * 2 - 1 dof(1) * 2 dof(2) * 2 - 1 dof(2) * 2];
T = [l m];
p(j) = E * A * T * U(index) / L;
end
end
- 空间桁架相关函数
function [S] = spacetruss(E, A, Wd, element, elementcon, xcord, ycord, zcord)
S = zeros(Wd, Wd);
for i = 1:element
dof = elementcon(i, :);
x = xcord(dof(2)) - xcord(dof(1));
y = ycord(dof(2)) - ycord(dof(1));
z = zcord(dof(2)) - zcord(dof(1));
L = sqrt(x^2 + y^2 + z^2);
l = x / L;
m = y / L;
n = z / L;
R = [l m n 0 0 0; 0 0 0 l m n];
k = (E * A / L) * [1 -1; -1 1];
K = R' * k * R;
index = [dof(1) * 3 - 2:dof(1) * 3 dof(2) * 3 - 2:dof(2) * 3];
S(index, index) = S(index, index) + K;
end
end
- `spacetrussforce`函数
function [p] = spacetrussforce(E, A, element, elementcon, xcord, ycord, zcord, U)
p = zeros(element, 1);
for j = 1:element
dof = elementcon(j, :);
x = xcord(dof(2)) - xcord(dof(1));
y = ycord(dof(2)) - ycord(dof(1));
z = zcord(dof(2)) - zcord(dof(1));
L = sqrt(x^2 + y^2 + z^2);
l = x / L;
m = y / L;
n = z / L;
index = [dof(1) * 3 - 2:dof(1) * 3 dof(2) * 3 - 2:dof(2) * 3];
T = [-l -m -n l m n];
p(j) = E * A * T * U(index) / L;
end
end
- 平面框架相关函数
function [S] = planeframe(E, A, I, Wd, element, elementcon, xcord, ycord)
S = zeros(Wd, Wd);
for j = 1:element
dof = elementcon(j, :);
x = xcord(dof(2)) - xcord(dof(1));
y = ycord(dof(2)) - ycord(dof(1));
L = sqrt(x^2 + y^2);
l = x / L;
m = y / L;
R = [l m 0 0 0 0; -m l 0 0 0 0; 0 0 1 0 0 0; 0 0 0 l m 0; 0 0 0 -m l 0; 0 0 0 0 0 1];
k = E * [A / L 0 0 -A / L 0 0; 0 12 * I / L^3 6 * I / L^2 0 -12 * I / L^3 6 * I / L^2; 0 6 * I / L^2 4 * I / L 0 -6 * I / L^2 2 * I / L; -A / L 0 0 A / L 0 0; 0 -12 * I / L^3 -6 * I / L^2 0 12 * I / L^3 -6 * I / L^2; 0 6 * I / L^2 2 * I / L 0 -6 * I / L^2 4 * I / L];
K = R' * k * R;
index = [3 * dof(1) - 2:3 * dof(1) 3 * dof(2) - 2:3 * dof(2)];
S(index, index) = S(index, index) + K;
end
end
- `drawplaneframe`函数
function drawplaneframe(element, elementcon, nodecord, newcord, newcord)
a = max(Snewcord(:, 1));
b = min(Snewcord(:, 1));
c = max(Snewcord(:, 2));
d = min(Snewcord(:, 2));
x0 = b * 1.3 - 1;
x1 = a * 1.3 + 1;
y0 = d * 1.3 - 1;
y1 = c * 1.3 + 1;
for n = 1:element
dof = elementcon(n, :);
xx = [nodecord(dof(1), 1), nodecord(dof(2), 1)];
yy = [nodecord(dof(1), 2), nodecord(dof(2
### 第2章：框架结构的矩阵分析（续）
1. **平面框架相关函数（续）**
```matlab
function drawplaneframe(element, elementcon, nodecord, newcord, newcord)
a = max(Snewcord(:, 1));
b = min(Snewcord(:, 1));
c = max(Snewcord(:, 2));
d = min(Snewcord(:, 2));
x0 = b * 1.3 - 1;
x1 = a * 1.3 + 1;
y0 = d * 1.3 - 1;
y1 = c * 1.3 + 1;
for n = 1:element
dof = elementcon(n, :);
xx = [nodecord(dof(1), 1), nodecord(dof(2), 1)];
yy = [nodecord(dof(1), 2), nodecord(dof(2), 2)];
xxn = [Snewcord(dof(1), 1), Snewcord(dof(2), 1)];
yyn = [Snewcord(dof(1), 2), Snewcord(dof(2), 2)];
plot(xx, yy, '-.k', xxn, yyn, 'k');
grid on;
title('Scaled Deformation');
xlabel('Horizontal Deformation');
ylabel('Vertical deformation');
axis([x0, x1, y0, y1]);
hold on;
end
end
- 空间框架相关函数
function [S] = spaceframe(E, A, Ge, IZ, IY, J, Wd, element, elementcon, xcord, ycord, zcord, beta)
S = zeros(Wd, Wd);
for j = 1:element
dof = elementcon(j, :);
x = xcord(dof(2)) - xcord(dof(1));
y = ycord(dof(2)) - ycord(dof(1));
z = zcord(dof(2)) - zcord(dof(1));
if xcord(dof(2)) == xcord(dof(1)) && ycord(dof(2)) == ycord(dof(1))
if zcord(dof(2)) > zcord(dof(1))
T = [0 0 1; 0 1 0; -1 0 0];
else
T = [0 0 -1; 0 1 0; 1 0 0];
end
else
L = sqrt(x^2 + y^2 + z^2);
l = x / L;
m = y / L;
n = z / L;
D = sqrt(l^2 + m^2);
T = [l m n; -m / D l / D 0; -l * n / D -m * n / D D];
end
b = beta(j);
bb = [1 0 0; 0 cos(b) sin(b); 0 -sin(b) cos(b)];
ZZ = zeros(3, 3);
B = [bb ZZ ZZ ZZ; ZZ bb ZZ ZZ; ZZ ZZ bb ZZ; ZZ ZZ ZZ bb];
R = [T ZZ ZZ ZZ; ZZ T ZZ ZZ; ZZ ZZ T ZZ; ZZ ZZ ZZ T];
k1 = E * [A / L 0 0 0 0 0 0 12 * IZ / L^3 0 0 0 6 * IZ / L^2 0 0 12 * IY / L^3 0 -6 * IY / L^2 0; 0 0 0 Ge * J / (L * E) 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 6 * IZ / L^2 0 0 0 4 * IZ / L 0 0 -6 * IY / L^2 0 0 0 2 * IY / L 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
k2 = E * [-A / L 0 0 0 0 0 0 0 -6 * IY / L^2 0 4 * IY / L 0 0 6 * IZ / L^2 0 0 0 4 * IZ / L; 0 0 0 -Ge * J / (L * E) 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 -6 * IZ / L^2 0 0 0 2 * IZ / L 0 0 6 * IY / L^2 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
k3 = k2';
k4 = E * [A / L 0 0 0 0 0 0 -6 * IZ / L^2 0 0 0 2 * IZ / L; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 6 * IZ / L^2 0 0 0 4 * IZ / L 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
k = [k1 k2; k3 k4];
K = R' * B' * k * B * R;
index = [6 * dof(1) - 5:6 * dof(1) 6 * dof(2) - 5:6 * dof(2)];
S(index, index) = S(index, index) + K;
end
end
- `drawspaceframe`函数
function drawspaceframe(element, elementcon, nodecord, newcord, Snewcord)
for n = 1:element
dof = elementcon(n, :);
xx = [nodecord(dof(1), 1), nodecord(dof(2), 1)];
yy = [nodecord(dof(1), 2), nodecord(dof(2), 2)];
zz = [nodecord(dof(1), 3), nodecord(dof(2), 3)];
xxn = [Snewcord(dof(1), 1), Snewcord(dof(2), 1)];
yyn = [Snewcord(dof(1), 2), Snewcord(dof(2), 2)];
zzn = [Snewcord(dof(1), 3), Snewcord(dof(2), 3)];
plot3(zz, xx, yy, '-.k', zzn, xxn, yyn, 'k');
grid on;
title('Scaled Deformation');
hold on;
end
end
- 网格相关函数
function [S] = grid(E, I, Ge, J, Wd, element, elementcon, xcord, ycord)
S = zeros(Wd, Wd);
for j = 1:element
dof = elementcon(j, :);
x = xcord(dof(2)) - xcord(dof(1));
y = ycord(dof(2)) - ycord(dof(1));
L = sqrt(x^2 + y^2);
l = x / L;
m = y / L;
R = [1 0 0 0 0 0; 0 -l -m 0 0 0; 0 m -l 0 0 0; 0 0 0 1 0 0; 0 0 0 0 -l -m; 0 0 0 0 m -l];
k1 = E * [12 * I / L^3 0 6 * I / L^2; 0 Ge * J / (E * L) 0; 6 * I / L^2 0 4 * I / L];
k2 = E * [-12 * I / L^3 0 6 * I / L^2; 0 -Ge * J / (E * L) 0; -6 * I / L^2 0 2 * I / L];
k3 = k2';
k4 = E * [12 * I / L^3 0 -6 * I / L^2; 0 Ge * J / (E * L) 0; -6 * I / L^2 0 4 * I / L];
k = [k1 k2; k3 k4];
K = R' * k * R;
index = [3 * dof(1) - 2:3 * dof(1) 3 * dof(2) - 2:3
### 第2章：框架结构的矩阵分析（续）
3. **网格相关函数（续）**
```matlab
function [S] = grid(E, I, Ge, J, Wd, element, elementcon, xcord, ycord)
S = zeros(Wd, Wd);
for j = 1:element
dof = elementcon(j, :);
x = xcord(dof(2)) - xcord(dof(1));
y = ycord(dof(2)) - ycord(dof(1));
L = sqrt(x^2 + y^2);
l = x / L;
m = y / L;
R = [1 0 0 0 0 0; 0 -l -m 0 0 0; 0 m -l 0 0 0; 0 0 0 1 0 0; 0 0 0 0 -l -m; 0 0 0 0 m -l];
k1 = E * [12 * I / L^3 0 6 * I / L^2; 0 Ge * J / (E * L) 0; 6 * I / L^2 0 4 * I / L];
k2 = E * [-12 * I / L^3 0 6 * I / L^2; 0 -Ge * J / (E * L) 0; -6 * I / L^2 0 2 * I / L];
k3 = k2';
k4 = E * [12 * I / L^3 0 -6 * I / L^2; 0 Ge * J / (E * L) 0; -6 * I / L^2 0 4 * I / L];
k = [k1 k2; k3 k4];
K = R' * k * R;
index = [3 * dof(1) - 2:3 * dof(1) 3 * dof(2) - 2:3 * dof(2)];
S(index, index) = S(index, index) + K;
end
end
- `drawgrid`函数
function drawgrid(element, elementcon, nodecord, Snewcord)
for n = 1:element
dof = elementcon(n, :);
xx = [nodecord(dof(1), 1), nodecord(dof(2), 1)];
zz = [nodecord(dof(1), 2), nodecord(dof(2), 2)];
yy = [nodecord(dof(1), 3), nodecord(dof(2), 3)];
xxn = [Snewcord(dof(1), 1), Snewcord(dof(2), 1)];
zzn = [Snewcord(dof(1), 2), Snewcord(dof(2), 2)];
yyn = [Snewcord(dof(1), 3), Snewcord(dof(2), 3)];
plot3(yy, zz, xx, 'k-.', yyn, zzn, xxn, 'k');
title('Scaled deformation');
xlabel('z');
ylabel('x');
zlabel('y');
hold on;
end
end
- 特殊情况处理相关函数
function [F] = nodeforces(Wd, eledist, loaddefiner, loaddist, xcord, ycord, option, F)
for k = 1:size(eledist, 1)
OPT = option(k, :);
dof = eledist(k, :);
x = xcord(dof(2)) - xcord(dof(1));
y = ycord(dof(2)) - ycord(dof(1));
L = sqrt(x^2 + y^2);
l = x / L;
m = y / L;
R = [l m 0 0 0 0; -m l 0 0 0 0; 0 0 1 0 0 0; 0 0 0 l m 0; 0 0 0 -m l 0; 0 0
### 第2章：框架结构的矩阵分析（续）
4. **特殊情况处理相关函数（续）**
```matlab
function [F] = nodeforces(Wd, eledist, loaddefiner, loaddist, xcord, ycord, option, F)
for k = 1:size(eledist, 1)
OPT = option(k, :);
dof = eledist(k, :);
x = xcord(dof(2)) - xcord(dof(1));
y = ycord(dof(2)) - ycord(dof(1));
L = sqrt(x^2 + y^2);
l = x / L;
m = y / L;
R = [l m 0 0 0 0; -m l 0 0 0 0; 0 0 1 0 0 0; 0 0 0 l m 0; 0 0 0 -m l 0; 0 0 0 0 0 1];
switch OPT
case 'dist'
a1 = loaddefiner(k, 1);
a2 = loaddefiner(k, 2);
b = loaddefiner(k, 3);
dist = loaddist(k, :)';
c1 = b * (L - a1 - b / 2) / L;
c2 = (2 * b * L - 2 * ((a1 + b)^3 - a1^3) / L + ((a1 + b)^4 - a1^4) / L^2) / (2 * L);
c3 = (6 * L^2 * ((a1 + b)^2 - a1^2) - 8 * L * ((a1 + b)^3 - a1^3) + 3 * ((a1 + b)^4 - a1^4)) / (12 * L^2);
c4 = b * (L - a2 - b / 2) / L;
c5 = (2 * b * L - 2 * ((a2 + b)^3 - a2^3) / L + ((a2 + b)^4 - a2^4) / L^2) / (2 * L);
c6 = (6 * L^2 * ((a2 + b)^2 - a2^2) - 8 * L * ((a2 + b)^3 - a2^3) + 3 * ((a2 + b)^4 - a2^4)) / (12 * L^2);
case 'conc'
a1 = loaddefiner(k, 1);
a2 = loaddefiner(k, 2);
b = loaddefiner(k, 3);
dist = loaddist(k, :)';
c1 = (L - a1) / L;
c2 = ((L - a1)^2 * (L + 2 * a1)) / L^3;
c3 = (a1 * (L - a1)^2) / L^2;
c4 = (L - a2) / L;
c5 = ((L - a2)^2 * (L + 2 * a2)) / L^3;
c6 = (a2 * (L - a2)^2) / L^2;
end
C = [-c1 0 0 0 0 0; 0 -c2 0 0 0 0; 0 -c3 0 0 0 0; 0 0 0 -c4 0 0; 0 0 0 0 -c5 0; 0 0 0 0 0 -c6];
pij = C * R * dist;
index = [3 * dof(1) - 2:3 * dof(1) 3 * dof(2) - 2:3 * dof(2)];
F(index) = F(index) + (-R' * pij);
end
end
- `settle`函数（处理支座沉降）
function [F] = settle(E, A, I, Wd, elementcon, settleelement, settlevalue, settlenode, xcord, ycord, F)
n = size(settleelement, 1);
for i = 1:n
### 第2章：框架结构的矩阵分析（续）
4. **特殊情况处理相关函数（续）**
```matlab
function [F] = settle(E, A, I, Wd, elementcon, settleelement, settlevalue, settlenode, xcord, ycord, F)
n = size(settleelement, 1);
for i = 1:n
connect = settleelement(i, :);
x = xcord(connect(2)) - xcord(connect(1));
y = ycord(connect(2)) - ycord(connect(1));
L = sqrt(x^2 + y^2);
l = x / L;
m = y / L;
kij = E * [-A / L 0 0; 0 -12 * I / L^3 6 * I / L^2; 0 -6 * I / L^2 2 * I / L];
R = [l m 0; -m l 0; 0 0 1];
if connect(1) == settlenode(i)
delta1 = settlevalue(i, :)';
p = kij' * R * delta1;
index = [3 * connect(2) - 2 3 * connect(2) - 1 3 * connect(2)];
else
delta2 = settlevalue(i, :)';
index = [3 * connect(1) - 2 3 * connect(1) - 1 3 * connect(1)];
p = kij * R * delta2;
end
P = -R' * p;
F(index) = F(index) + P;
end
end
- `temperature`函数（处理温度变化）
function [F] = temperature(E, A, I, h, alfa, Wd, tempelement, temperaturedeg, elementcon, xcord, ycord, F)
n = size(tempelement, 1);
for i = 1:n
T1 = temperaturedeg(i, 1);
T2 = temperaturedeg(i, 2);
connect = tempelement(i, :);
x = xcord(connect(1)) - xcord(connect(2));
y = ycord(connect(1)) - ycord(connect(2));
L = sqrt(x^2 + y^2);
l = x / L;
m = y / L;
Rij = [l m 0; -m l 0; 0 0 1];
Rji = [-l -m 0; m -l 0; 0 0 1];
pij = alfa * [-E * A * (T2 + T1) / 2; 0; E * I * (T1 - T2) / h];
pji = alfa * [-E * A * (T2 + T1) / 2; 0; E * I * (T2 - T1) / h];
Pi = -Rij' * pij;
Pj = -Rji' * pji;
P = [Pi; Pj];
index = [3 * connect(1) - 2:3 * connect(1) 3 * connect(2) - 2:3 * connect(2)];
F(index) = F(index) + P;
end
end
- `aberration`函数（处理构件不匹配）
function [F] = aberration(E, A, I, Wd, abelement, abvalue, xcord, ycord, option, F)
n = size(abelement, 1);
for i = 1:n
dof = abelement(i, :);
switch option
case 'truss'
index1 = [dof(1) * 2 - 1 dof(1) * 2];
index2 = [dof(2) * 2 - 1 dof(2) * 2];
x = xcord(dof(2)) - xcord(dof(1));
y = ycord(dof(2)) - ycord(dof(1));
L = sqrt(x^2 + y^2);
l = x / L;
m = y / L;
R = [l m];
k
### 第2章：框架结构的矩阵分析（续）
4. **特殊情况处理相关函数（续）**
```matlab
function [F] = aberration(E, A, I, Wd, abelement, abvalue, xcord, ycord, option, F)
n = size(abelement, 1);
for i = 1:n
dof = abelement(i, :);
switch option
case 'truss'
index1 = [dof(1) * 2 - 1 dof(1) * 2];
index2 = [dof(2) * 2 - 1 dof(2) * 2];
x = xcord(dof(2)) - xcord(dof(1));
y = ycord(dof(2)) - ycord(dof(1));
L = sqrt(x^2 + y^2);
l = x / L;
m = y / L;
R = [l m];
kii = (E * A / L);
kij = -kii;
delta = abvalue(i);
case 'frame2d'
index1 = [3 * dof(1) - 2 3 * dof(1) - 1 3 * dof(1)];
index2 = [3 * dof(2) - 2 3 * dof(2) - 1 3 * dof(2)];
x = xcord(dof(2)) - xcord(dof(1));
y = ycord(dof(2)) - ycord(dof(1));
L = sqrt(x^2 + y^2);
l = x / L;
m = y / L;
R = [l m 0; -m l 0; 0 0 1];
kii = E * [A / L 0 0; 0 12 * I / L^3 6 * I / L^2; 0 6 * I / L^2 4 * I / L];
kij = E * [-A / L 0 0; 0 -12 * I / L^3 6 * I / L^2; 0 -6 * I / L^2 2 * I / L];
delta = [abvalue(i) 0 0];
end
pij = kii * delta;
pji = kij * delta;
Pi = -R' * pij;
Pj = -R' * pji;
F(index1) = F(index1) + Pi;
F(index2) = F(index2) + Pj;
end
end
- `discon`函数（处理释放端条件）
function [S] = discon(E, A, I, Wd, element, elementcon, xcord, ycord, disc, option)
S = zeros(Wd, Wd);
n = size(disc, 1);
for j = 1:element - n
dof = elementcon(j, :);
index = [3 * dof(1) - :3 * dof(1) 3 * dof(2) - :3 * dof(2)];
x = xcord(dof(2)) - xcord(dof(1));
y = ycord(dof(2)) - ycord(dof(1));
L = sqrt(x^2 + y^2);
l = x / L;
m = y / L;
R = [l m 0 0 0 0; -m l 0 0 0 0; 0 0 1 0 0 0; 0 0 0 l m 0; 0 0 0 -m l 0; 0 0 0 0 0 1];
k = E * [A / L 0 0 -A / L 0 0; 0 12 * I / L^3 6 * I / L^2 0 -12 * I / L^3 6 * I / L^2; 0 6 * I / L^2 4 * I / L 0 -6 * I / L^2 2 * I / L; -A / L 0 0 A / L 0 0; 0 -12 *
### 第2章：框架结构的矩阵分析（续）
4. **特殊情况处理相关函数（续）**
```matlab
function [S] = discon(E, A, I, Wd, element, elementcon, xcord, ycord, disc, option)
S = zeros(Wd, Wd);
n = size(disc, 1);
for j = 1:element - n
dof = elementcon(j, :);
index = [3 * dof(1) - :3 * dof(1) 3 * dof(2) - :3 * dof(2)];
x = xcord(dof(2)) - xcord(dof(1));
y = ycord(dof(2)) - ycord(dof(1));
L = sqrt(x^2 + y^2);
l = x / L;
m = y / L;
R = [l m 0 0 0 0; -m l 0 0 0 0; 0 0 1 0 0 0; 0 0 0 l m 0; 0 0 0 -m l 0; 0 0 0 0 0 1];
k = E * [A / L 0 0 -A / L 0 0; 0 12 * I / L^3 6 * I / L^2 0 -12 * I / L^3 6 * I / L^2; 0 6 * I / L^2 4 * I / L 0 -6 * I / L^2 2 * I / L; -A / L 0 0 A / L 0 0; 0 -12 * I / L^3 -6 * I / L^2 0 12 * I / L^3 -6 * I / L^2; 0 6 * I / L^2 2 * I / L 0 -6 * I / L^2 4 * I / L];
K = R' * k * R;
S(index, index) = S(index, index) + K;
end
for ii = 1:n
dof = disc(ii, :);
OPT = option(ii, :);
index = [3 * dof(1) - 2 3 * dof(1) - 1 3 * dof(1) 3 * dof(2) - 2 3 * dof(2) - 1 3 * dof(2)];
x = xcord(dof(2)) - xcord(dof(1));
y = ycord(dof(2)) - ycord(dof(1));
L = sqrt(x^2 + y^2);
l = x / L;
m = y / L;
R = [l m 0 0 0 0; -m l 0 0 0 0; 0 0 1 0 0 0; 0 0 0 l m 0; 0 0 0 -m l 0; 0 0 0 0 0 1];
switch OPT
case 'pin'
k1 = E * [A / L 0 0; 0 3 * I / L^3 0; 0 0 0];
k2 = E * [A / L 0 0; 0 -3 * I / L^3 3 * I / L^2; 0 0 0];
k3 = E * [A / L 0 0; 0 3 * I / L^3 -3 * I / L^2; 0 -3 * I / L^2 3 * I / L];
k = [k1 k2; k2' k3];
case 'trs'
k1 = E * [A / L 0 0; 0 0 0; 0 0 0];
k2 = E * [-A / L 0 0; 0 0 0; 0 0 0];
k3 = E * [A / L 0 0; 0 0 0; 0 0 0];
k = [k1 k2; k2' k3];
case 'shr'
k1 = E * [A
### 第3章：使用有限元程序进行结构的弹性分析
1. **获取高斯积分点坐标和权重的函数**
```matlab
function [point, weight] = gausslegendre(number)
switch number
case 'one'
point = [0];
weight = [2];
case 'two'
point = [-0.577350269189626; 0.577350269189626];
weight = [1; 1];
case 'three'
point = [-.774596669241483; 0;.774596669241483];
weight = [.555555555555555;.888888888888888;.555555555555555];
case 'four'
point = [-.861136311594053; -.339981043584856;.339981043584856;.861136311594053];
weight = [.347854845137454;.652145154862546;.652145154862546;.347854845137454];
end
end
- 计算一维单元插值（形状）函数及其导数的函数
function [shape, RDerivatives] = shapeR(r, option)
switch option
case 'twonode'
shape = [(1 - r) (1 + r)] / 2;
RDerivatives = [-1 1] / 2;
case 'threenode'
shape = [(1 - r) - (1 - r^2) (1 + r) - (1 - r^2) 2 * (1 - r^2)] / 2;
RDerivatives = [2 * r - 1 2 * r + 1 -4 * r] / 2;
case 'fournode'
shape = [(1 - r) - (1 - r^2) + (-9 * r^3 + r^2 + 9 * r - 1) / 8 (1 + r) - (1 - r^2) + (9 * r^3 + r^2 - 9 * r - 1) / 8 2 * (1 - r^2) + (27 * r^3 + 7 * r^2 - 27 * r - 7) / 8 (-27 * r^3 - 9 * r^2 + 27 * r + 9) / 8] / 2;
RDerivatives = [(9 * r) / 4 - (27 * r^2) / 8 + 1 / 8 (27 * r^2) / 8 + (9 * r) / 4 - 1 / 8 (81 * r^2) / 8 - (9 * r) / 4 - 27 / 8 27 / 8 - (81 * r^2) / 8 - (9 * r) / 4] / 2;
end
end
- 计算二维结构刚度矩阵（平面应力、平面应变）的函数
function [S] = continuum2d(Wd, C, N, NodeNumbertotal, thickness, nodes, cord, gausepoint, weight, option)
S = zeros(Wd, Wd);
for i = 1:N
dof = nodes(i, :);
eledof = [dof dof + NodeNumbertotal];
betaLength = length(dof);
for j = 1:size(gausepoint, 1)
point = gausepoint(j, :);
r = point(1);
s = point(2);
[shape, RSderivation] = shapeStress(r, s, option);
elecord = cord(dof, :);
[XYderivation, Jacob
### 第4章：使用有限元程序进行结构的弹塑性分析
1. **计算变形梯度、拉格朗日应变和第二Piola - Kirchhoff应力的相关代码片段（示例）**
- 计算变形梯度（以平面应变条件下的弯曲和牵引荷载为例，在主程序中可能的计算方式）
```matlab
% 假设已经定义了节点坐标node和位移向量U
for i = 1:length(node)
x(i) = node(i, 1) + U(2 * i - 1);
y(i) = node(i, 2) + U(2 * i);
end
% 计算变形梯度F
for i = 1:length(element)
node1 = element(i, 1);
node2 = element(i, 2);
dx_dX = (x(node2) - x(node1)) / (node(node2, 1) - node(node1, 1));
dy_dX = (y(node2) - y(node1)) / (node(node2, 1) - node(node1, 1));
dx_dY = 0; % 平面应变假设
dy_dY = 1; % 平面应变假设
F(i, :, :) = [dx_dX dx_dY; dy_dX dy_dY];
end
E = zeros(length(element), 3);
for i = 1:length(element)
F_i = F(i, :, :);
E(i, 1) = 0.5 * (F_i(1, 1)^2 + F_i(2, 1)^2 - 1); % E_xx
E(i, 2) = 0.5 * (F_i(1, 2)^2 + F_i(2, 2)^2 - 1); % E_yy
E(i, 3) = F_i(1, 1) * F_i(1, 2) + F_i(2, 1) * F_i(2, 2); % E_xy
end
- 计算第二Piola - Kirchhoff应力（假设材料参数已定义，如弹性模量E和泊松比nu，以平面应力为例）
C = E / (1 - nu^2) * [1 nu 0; nu 1 0; 0 0 (1 - nu) / 2]; % 平面应力本构矩阵
S = zeros(length(element), 3);
for i = 1:length(element)
E_i = [E(i, 1); E(i, 2); E(i, 3)];
S(i, :) = C * E_i;
end
- 应力积分相关函数（以结合等向/随动硬化的材料为例，可能涉及的计算过程，未完整封装函数）
- 计算等效塑性应变增量（根据公式计算，可能在迭代过程中调用）
function [DeltaEpsilonP] = calculateDeltaEpsilonP(hatSigmaPr, SigmaY, G, C, M)
% 计算混合硬化斜率
hatC = (1 - M) * C;
% 计算函数f(DeltaEpsilonP)的值
f = @(DeltaEpsilonP) hatSigmaPr / (SigmaY + (3 / 2) * (2 * G + hatC) * DeltaEpsilonP) - 1;
% 使用牛顿 - 拉夫森法迭代求解DeltaEpsilonP
DeltaEpsilonP = newtonRaphsonSolve(f, initialGuess, tolerance);
end
- 更新应力分量（根据计算得到的等效塑性应变增量等参数更新应力）
function [sigma] = updateStress(hatSigmaPr, SigmaY, eta, alpha0, DeltaAlpha, G, C, M)
hatC = (1 - M) * C;
% 根据公式计算新的应力分量
sigma = zeros(6, 1); % 假设应力向量维度为6（根据具体问题调整）
for i = 1:6
sigma(i) = sigmaPr(i) + eta(i) * SigmaY + alpha0(i) + DeltaAlpha(i);
end
end
- 求解非线性问题相关代码片段（如牛顿 - 拉夫森迭代法的应用示例，在主程序中可能的实现方式）
第5章：有限变形与超弹性
- 计算变形梯度、拉格朗日应变和第二Piola - Kirchhoff应力的相关代码片段（示例）
- 计算变形梯度（以平面应变条件下的弯曲和牵引荷载为例，在主程序中可能的计算方式）
% 假设已经定义了节点坐标node和位移向量U
for i = 1:length(node)
x(i) = node(i, 1) + U(2 * i - 1);
y(i) = node(i, 2) + U(2 * i);
end
% 计算变形梯度F
for i = 1:length(element)
node1 = element(i, 1);
node2 = element(i, 2);
dx_dX = (x(node2) - x(node1)) / (node(node2, 1) - node(node1, 1));
dy_dX = (y(node2) - y(node1)) / (node(node2, 1) - node(node1, 1));
dx_dY = 0; % 平面应变假设
dy_dY = 1; % 平面应变假设
F(i, :, :) = [dx_dX dx_dY; dy_dX dy_dY];
end
E = zeros(length(element), 3);
for i = 1:length(element)
F_i = F(i, :, :);
E(i, 1) = 0.5 * (F_i(1, 1)^2 + F_i(2, 1)^2 - 1); % E_xx
E(i, 2) = 0.5 * (F_i(1, 2)^2 + F_i(2, 2)^2 - 1); % E_yy
E(i, 3) = F_i(1, 1) * F_i(1, 2) + F_i(2, 1) * F_i(2, 2); % E_xy
end
- 计算第二Piola - Kirchhoff应力（假设材料参数已定义，如弹性模量E和泊松比nu，以平面应力为例）
C = E / (1 - nu^2) * [1 nu 0; nu 1 0; 0 0 (1 - nu) / 2]; % 平面应力本构矩阵
S = zeros(length(element), 3);
for i = 1:length(element)
E_i = [E(i, 1); E(i, 2); E(i, 3)];
S(i, :) = C * E_i;
end
- 超弹性材料拟合相关代码片段（以neo - Hookean模型为例，可能的实现方式）
function [stress] = neoHookeanModel(lambda, mu, F)
% 计算变形梯度的行列式
J = det(F);
% 计算右柯西 - 格林变形张量
C = F' * F;
% 计算不变量
I1 = trace(C);
% 计算应力
stress = 2 * mu * (F - 1 / 2 * I1 * eye(3)) + lambda * log(J) * eye(3);
end
- 进行数据拟合（假设实验数据为stressData和stretchData）
% 初始参数猜测
lambda0 = 100; % 初始拉梅常数lambda猜测值
mu0 = 50; % 初始剪切模量mu猜测值
% 定义优化选项
options = optimset('MaxIter', 1000, 'TolFun', 1e-6);
% 使用最小二乘法拟合模型参数
[params, residual] = lsqcurvefit(@neoHookeanModel, [lambda0, mu0], stretchData, stressData, [], [], options);
lambda_fitted = params(1);
mu_fitted = params(2);
第6章：有限应变
- 基于特定理论（如TL和UL公式）的有限应变分析相关代码片段（示例）
- 计算变形梯度（以平面应变条件下的弯曲和牵引荷载为例，在主程序中可能的计算方式，与第5章类似但可能涉及不同的处理或计算顺序）
第6章：有限应变（续）
- 基于特定理论（如TL和UL公式）的有限应变分析相关代码片段（示例）
- 计算变形梯度（以平面应变条件下的弯曲和牵引荷载为例，在主程序中可能的计算方式，与第5章类似但可能涉及不同的处理或计算顺序）
% 假设已经定义了节点坐标node和位移向量U
for i = 1:length(node)
x(i) = node(i, 1) + U(2 * i - 1);
y(i) = node(i, 2) + U(2 * i);
end
% 计算变形梯度F
for i = 1:length(element)
node1 = element(i, 1);
node2 = element(i, 2);
dx_dX = (x(node2) - x(node1)) / (node(node2, 1) - node(node1, 1));
dy_dX = (y(node2) - y(node1)) / (node(node2, 1) - node(node1, 1));
dx_dY = 0; % 平面应变假设
dy_dY = 1; % 平面应变假设
F(i, :, :) = [dx_dX dx_dY; dy_dX dy_dY];
end
- 计算拉格朗日应变（基于TL公式，使用上述计算的变形梯度F）
E_TL = zeros(length(element), 3);
for i = 1:length(element)
F_i = F(i, :, :);
E_TL(i, 1) = 0.5 * (F_i(1, 1)^2 + F_i(2, 1)^2 - 1); % E_xx
E_TL(i, 2) = 0.5 * (F_i(1, 2)^2 + F_i(2, 2)^2 - 1); % E_yy
E_TL(i, 3) = F_i(1, 1) * F_i(1, 2) + F_i(2, 1) * F_i(2, 2); % E_xy
end
- 计算拉格朗日应力（假设材料参数已定义，如弹性模量E和泊松比nu，以平面应力为例，使用TL公式计算变形梯度F）
C = E / (1 - nu^2) * [1 nu 0; nu 1 0; 0 0 (1 - nu) / 2]; % 平面应力本构矩阵
S_TL = zeros(length(element), 3);
for i = 1:length(element)
E_i = [E_TL(i, 1); E_TL(i, 2); E_TL(i, 3)];
S_TL(i, :) = C * E_i;
end
- 计算更新的拉格朗日（UL）公式相关量（如UL应变和应力，假设在计算变形梯度F后进行）
% 计算UL变形梯度F_UL（假设已经有了F）
F_UL = zeros(size(F));
for i = 1:length(element)
F_UL(i, :, :) = inv(F(i, :, :))' * F(i, :, :);
end
% 计算UL应变E_UL
E_UL = zeros(length(element), 3);
for i = 1:length(element)
F_UL_i = F_UL(i, :, :);
E_UL(i, 1) = 0.5 * (F_UL_i(1, 1)^2 + F_UL_i(2, 1)^2 - 1); % E_xx
E_UL(i, 2) = 0.5 * (F_UL_i(1, 2)^2 + F_UL_i(2, 2)^2 - 1); % E_yy
E_UL(i, 3) = F_UL_i(1, 1) * F_UL_i(1, 2) + F_UL_i(2, 1) * F_UL_i(2, 2); % E_xy
end
% 计算UL应力S_UL（假设材料参数不变）
S_UL = zeros(length(element), 3);
for i = 1:length(element)
E_UL_i = [E_UL(i,
### 第7章：线性方程组和特征问题的求解
1. **线性方程组求解相关函数**
- `mldivide`函数：这是MATLAB中用于求解线性方程组`Ax = B`的内置函数，其用法为`x = A\B`，其中`A`是系数矩阵，`B`是右端项向量，`x`是解向量。例如，对于方程组`A = [1 2; 3 4]; B = [5; 6]; x = A\B`，将计算出满足`Ax = B`的`x`值。
- `qr`函数：用于对矩阵进行QR分解，语法为`[Q, R] = qr(A)`，其中`A`是待分解的矩阵，`Q`是正交矩阵，`R`是上三角矩阵。例如，`A = [1 2 3; 4 5 6; 7 8 9]; [Q, R] = qr(A)`将对矩阵`A`进行QR分解并返回`Q`和`R`。
- `lu`函数：执行LU分解，语法为`[L, U, P] = lu(A)`，`A`是输入矩阵，`L`是下三角矩阵（主对角线元素为1），`U`是上三角矩阵，`P`是置换矩阵。例如，`A = [2 1; 4 3]; [L, U, P] = lu(A)`将得到矩阵`A`的LU分解结果。
- `chol`函数：进行Cholesky分解，语法为`R = chol(A)`，其中`A`必须是对称正定矩阵，`R`是上三角矩阵，满足`A = R'*R`。例如，`A = [4 2; 2 2]; R = chol(A)`将计算矩阵`A`的Cholesky分解。
2. **特征问题求解相关函数**
- `eig`函数：用于计算矩阵的特征值和特征向量，语法为`[V, D] = eig(A)`，其中`A`是输入矩阵，`V`是特征向量矩阵，`D`是对角矩阵，其对角元素是特征值。例如，`A = [1 2; 3 4]; [V, D] = eig(A)`将返回矩阵`A`的特征值和特征向量。
- `jacobi`方法相关函数（假设自行编写的`jacobi`函数实现雅可比方法）
```matlab
function [V, D] = jacobi(A, tol)
% 初始化特征向量矩阵为单位矩阵
V = eye(size(A));
% 计算矩阵的对角元素
D = diag(diag(A));
% 计算非对角元素的平方和
off = sum(sum(A.^2)) - sum(diag(A).^2);
% 迭代直到收敛
while off > tol
% 找到最大非对角元素的位置
[p, q] = find(abs(A - D) == max(max(abs(A - D))));
p = p(1);
q = q(1);
% 计算旋转角度
if A(p, q) ~= 0
theta = 0.5 * atan(2 * A(p, q) / (D(q, q) - D(p, p)));
else
theta = 0;
end
% 计算旋转矩阵
c = cos(theta);
s = sin(theta);
J = [c s; -s c];
% 更新矩阵A和特征向量矩阵V
A([p, q], :) = J' * A([p, q], :);
A(:, [p, q]) = A(:, [p, q]) * J;
V(:, [p, q]) = V(:, [p, q]) * J;
% 更新对角矩阵D
D([p, q], [p, q]) = J' * D([p, q], [p, q]) * J;
% 更新非对角元素的平方和
off = sum(sum(A.^2)) - sum(diag(A).^2);
end
end
generalized jacobi方法相关函数（假设自行编写的generalized_jacobi函数实现广义雅可比方法，此处仅为示意，实际函数可能更复杂）
function [V, D] = generalized_jacobi(A, B, tol
浙公网安备 33010602011771号