MATLAB基础运算【5】：图像绘制学科专题

系列更新MATLAB各种运算和操作的说明。

本篇主要包括图像绘制+学科专题相关内容（可利用侧边目录直接定位）：

1. 非稳态扩散现象浓度分布函数
2. 管道层流的雷诺数与水头损失系数关系近似
3. 振动曲线方程
4. 强迫振动在简谐激励作用下对系统响应振幅影响
5. 半无限大物体的导热问题中的温度分布
6. 温度分布函数

前文：

MATLAB基础运算【1】：矩阵元素

MATLAB基础运算【2】：矩阵运算

MATLAB基础运算【3】：微积分运算

MATLAB基础运算【4】：图像绘制

1- 非稳态扩散现象浓度分布函数

已知非稳态扩散现象浓度分布函数，分别使用笛卡尔直角坐标及半对数坐标绘制下列两种情况下的函数曲线图（0≤x≤20）：

1）D=1, t分别取1, 2, 4, 8, 16时u的分布图（绘制在同一张图上）;

2）D=100, t分别取1, 2, 4, 8, 16, 32时u的分布图。

D=1，半对数，命令行输入：

x = 0:0.05:20; % 创建x的向量，从0开始，以0.05为步长，直到20
t = [1, 2, 4, 8, 16];
for i = 1:length(t) % 循环遍历t中的每个时间值
u = (1 ./ sqrt(4. * pi .* t(i))) .* exp(-x .^ 2 ./ (4. * t(i)));
semilogy(x, u); % 绘制u的对数-对数图，使用semilogy函数
hold on % 保持当前图形，以便在同一个坐标轴上绘制下一条曲线
end
xlabel('x');
ylabel('u(semilogy)');
title('D=1');
legend('t=1','t=2','t=4','t=8','t=16') % 图例

输出图像：

 D=1，笛卡尔，命令行输入：

x = 0:0.05:20; % 创建x的向量，从0开始，以0.05为步长，直到20
t = [1, 2, 4, 8, 16];
for i = 1:length(t) % 循环遍历t中的每个时间值
u = (1 ./ sqrt(4. * pi .* t(i))) .* exp(-x .^ 2 ./ (4. * t(i)));
plot(x, u);
hold on % 保持当前图形，以便在同一个坐标轴上绘制下一条曲线
end
xlabel('x');
ylabel('u(linear)');
title('D=1');
legend('t=1','t=2','t=4','t=8','t=16') % 图例

输出图像：

 D=100，半对数，命令行输入：

x = 0:0.05:20; % 创建x的向量，从0开始，以0.05为步长，直到20
t = [1, 2, 4, 8, 16];
for i = 1:length(t) % 循环遍历t中的每个时间值
u = (1 ./ sqrt(4. * pi .* t(i).*100)) .* exp(-x .^ 2 ./ (4.*100.*t(i)));
semilogy(x, u); % 绘制u的对数-对数图，使用semilogy函数
hold on % 保持当前图形，以便在同一个坐标轴上绘制下一条曲线
end
xlabel('x');
ylabel('u(semilogy)');
title('D=100');
legend('t=1','t=2','t=4','t=8','t=16') % 图例

输出图像：

  D=100，笛卡尔，命令行输入：

x = 0:0.05:20; % 创建x的向量，从0开始，以0.05为步长，直到20
t = [1, 2, 4, 8, 16];
for i = 1:length(t) % 循环遍历t中的每个时间值
u = (1 ./ sqrt(4. * pi .*100.*t(i))) .* exp(-x .^ 2 ./ (4.*100.* t(i)));
plot(x, u);
hold on % 保持当前图形，以便在同一个坐标轴上绘制下一条曲线
end
xlabel('x');
ylabel('u(linear)');
title('D=100');
legend('t=1','t=2','t=4','t=8','t=16') % 图例

输出图像：

2- 管道层流的雷诺数与水头损失系数关系近似

下表为某管道局部结构中层流流动的雷诺数Re与水头损失系数K的数据，一般可通过下列两式来近似K=K(Re)这一变化关系：

（1） ；（一阶近似）

（2）（高阶近似）（可以采用fit或lsqcurvefit，只需保证所用函数和结果一致即可）

分别用上述两式拟合得到未知系数并作出曲线（数据点用红色空心圆点，一阶拟合曲线用绿色实线，高阶拟合曲线用黑色点划线）。

Re	2	4	8	16	32	64	128	256	512
K	56.3615	28.2406	14.2390	7.34264	4.0626	2.6171	2.0271	1.7572	1.5907

命令行输入：

Re = [2, 4, 8, 16, 32, 64, 128, 256, 512];
K = [56.3615,28.2406,14.2390,7.34264,4.0626,2.6171,2.0271,1.7572,1.5907];
g = @(y, Re) (y(1) + y(2) ./ Re); % 定义函数g，它是一个匿名函数，接受y和Re作为输入
f = @(x, Re) ((x(1).^20) + (x(2) ./ Re).^20).^(1 ./ 20); % 定义函数f，它是一个匿名函数，接受x和Re作为输入
x = [10, 10];
y = [10, 10]; % 定义初始猜测值x和y
b = lsqcurvefit(g, y, Re, K); % 使用lsqcurvefit函数拟合g函数，以找到最佳拟合参数b
a = lsqcurvefit(f, x, Re, K); % 使用lsqcurvefit函数拟合f函数，以找到最佳拟合参数a
scatter(Re, K, 'r'); % 绘制K关于Re的散点图
x1 = 2:0.1:512;
x2 = 2:10:512;
hold on
z1 = g(b, x1);% 计算g函数在b参数和x1上的值
plot(x1, z1, '-g'); % 绘制z1关于x1的曲线，使用绿色实线
z2 = f(a, x2); % 计算f函数在a参数和x2上的值
plot(x2, z2, '.k'); % 绘制z2关于x2的点图，使用黑色点
xlabel('Re');
ylabel('K');
hold off;

输出图像：

3- 振动曲线方程

设振动曲线方程为

1）利用plot函数画出方程在（0,12）区间内的曲线，步长取0.01，横坐标为＇time＇，纵坐标为＇displacement＇。

2）用max和min函数求出曲线的最大值和最小值及其位置（t_max,y_max）,(t_min,y_min)，在曲线上用红色，大小12的圆点标出(‘r.’,’MarkerSize’,12)，并在(t_max+0.2，y_max+0.02)和(t_max+0.2，y_max-0.02)的位置上用text函数标出最大最小值点的坐标。

3）设为符号函数，用diff函数求出S的导数S1，并用fplot画出导数的曲线，观察零点位置。

4)利用int函数积分S1，积分限为（t_max,t_min）,与(2)中y_min-y_max的结果比较。

1）命令行输入：

t = 0:0.01:12;
y = exp(-3. * t) .* sin(4. * t);
plot(t, y);
xlabel('time');
ylabel('displacement');

输出图像：

2）命令行输入：

[y_max, t_max1] = max(y);
[y_min, t_min1] = min(y);
t_max = 0.01 * (t_max1 - 1);
t_min = (t_min1 - 1) * 0.01;
hold on
plot(t_min, y_min, 'r.', 'Markersize', 12); % 在t_min处绘制y_min的点，并设置标记大小
plot(t_max, y_max, 'r.', 'Markersize', 12); % 在t_max处绘制y_max的点，并设置标记大小
text(t_max + 0.2, y_max + 0.02, '(0.23,0.3991)'); % 在t_max的右侧添加文本，显示最大值点的坐标
text(t_min + 0.2, y_min - 0.02, '(1.02,-0.0378)'); % 在t_min的左侧添加文本，显示最小值点的坐标
hold off

输出图像：

3）命令行输入：

syms s t
s = exp(-3. * t) .* sin(4. * t);
s1 = diff(s, t, 1);
fplot(s1, [0, 10]);
xlabel('t');
ylabel('sl');
eqn = s1 == 0;
[tq, param, cond] = solve(eqn, t, 'ReturnConditions', true); % 解方程s1=0，返回条件

输出图像：

[0,10]之间零点一共有13个，其中在0.23和1.02处(0后的第一个和第二个零点)取得该区间内的最大值和最小值。可以看到，随着t增大，函数是迅速衰减的。

4）命令行输入：

int(s1,t,t_min,t_max), y_max-y_min

得到结果：

对s1在（t_max,t_min）积分得到的结果与y_min-y_max数值相同。

4- 强迫振动在简谐激励作用下对系统响应振幅影响

，讨论不同 ζ 下的关系图

命令行输入：

lambda = 0:0.05:3;
e = [0, 0.125, 0.15, 0.2, 0.25, 0.5, 0.707, 1, 2];
for i = 1:length(e) % 开始一个循环，遍历e中的每个值
b = 1 ./ sqrt((1 - lambda.^2).^2 + 4 .* (e(i) .^ 2) .* (lambda.^2)); % 根据lamda和当前的e值计算b
plot(lambda, b);
hold on
end
ylim([0, 10]);
legend('e=0', 'e=0.125', 'e=0.15', 'e=0.2', 'e=0.25', 'e=0.5', 'e=0.707', 'e=1', 'e=2'); % 添加图例，对应不同的e值
xlabel('lambda');
ylabel('b');

输出图像：

解释：振动理论中，讨论强迫振动在简谐激励作用下对系统响应振幅影响的公式中，放大因子beta (用b代替)与频率比lambda、阻尼比ζ（matlab绘图中用e代替）有关，其公式如下图所示。频率比lambda为激励频率与系统固有频率之比，阻尼比ζ (e)为系统阻尼与临界阻尼之比。在阻尼比取不同取值时，放大因子随频率比的变化呈现不同模式，lambda接近1时为共振区，远大于1时为高频区，远小于1时为低频区。

5- 半无限大物体的导热问题中的温度分布

已知一端温度固定的半无限大物体的导热问题解析解为

其中，T为温度，x为到固定温度端的距离，t为时间，T0为初始温度，Tw为固定端温度，a为热扩散系数。取T0=0°C，Tw=30°C，a=10^-6 m2/s，x在[0,1]m之间，

1）在一张图上作出时间t=0s, 2000s, 5000s, 10000s, 20000s时的温度分布图，并标注图例.

命令行输入：

x = 0:0.01:1;
t = [0, 2000, 5000, 10000, 20000];
for i = 1:length(t)
T = 30 - 30 * erf(x ./ (2 * sqrt(10^(-6) * t(i))));
plot(x, T);
hold on
end
legend('t=0', 't=2000', 't=5000', 't=10000', 't=20000');
xlabel('x');
ylabel('T');

输出图像：

6- 温度函数分布

长为4，宽为2的长方形区域的中心（2,1）处的温度为100°C，区域内的温度分布满足：

将区域划分成100×50的网格，用surface和contour函数做出温度分布。

命令行输入：

x = linspace(0,4,100);
y = linspace(0,2,50);
[X,Y] = meshgrid(x,y); % 使用meshgrid函数生成网格矩阵X和Y，用于等高线图的坐标
T = 100.*exp(-3.*(X-2).^2 - (Y-1).^2);
contour(X,Y,T) % 使用contour函数绘制T的等高线图
surface(X,Y,T) % 使用contour函数绘制T的三维表面云图

得到图像：

【本篇完】