2024年数学建模A题(1——3问代码)
第一问 and 第二问
clc;
clear;
close all;
A = (0.55/(2*pi))^2;
R = @(theta) 0.55/(2*pi)*theta;
numBody = 223;
head_theta0 = 32*pi;
syms theta
T = theta*sqrt(theta^2+1) + log(sqrt(theta^2+1) + theta);
C = subs(T, theta, head_theta0);
front_offset = 0.275;
side_offset = 0.15;
theta_vals = linspace(0, head_theta0, 2000);
spiral_X = R(theta_vals).*cos(theta_vals);
spiral_Y = R(theta_vals).*sin(theta_vals);
figure;
axis equal;
grid on;
xlabel('X'); ylabel('Y');
title('龙身体随时间运动的动画');
collision_detected = false;
for t = 0:50:500
clf; hold on; axis equal; grid on;
xlabel('X'); ylabel('Y');
title(['时间 t = ', num2str(t, '%.7f'), ' 秒']);
plot(spiral_X, spiral_Y, 'k--', 'LineWidth', 1);
question = T == C - 4*pi/0.55*t;
head_theta = double(vpasolve(question, theta, [0, head_theta0]));
head_X = double(R(head_theta) * cos(head_theta));
head_Y = double(R(head_theta) * sin(head_theta));
center_X_list = head_X;
center_Y_list = head_Y;
X = head_theta;
plot(head_X, head_Y, 'ro', 'MarkerFaceColor','r');
bodyRects = {};
for i = 1:numBody
if i == 1
L = 2.86; rect_length = 3.41;
else
L = 1.65; rect_length = 2.2;
end
syms Y
eqn = A*(Y^2 + X^2) - L^2 == 2*A*Y*X*cos(Y - X);
solution = double(vpasolve(eqn, Y, [X, X+pi/2]));
L_X = R(solution) * cos(solution);
L_Y = R(solution) * sin(solution);
plot(L_X, L_Y, 'bo', 'MarkerFaceColor','b');
center_X_list = [center_X_list, L_X];
center_Y_list = [center_Y_list, L_Y];
dx = L_X - center_X_list(end-1);
dy = L_Y - center_Y_list(end-1);
len = sqrt(dx^2 + dy^2);
ux = dx / len; uy = dy / len;
vx = -uy; vy = ux;
A_point = [center_X_list(end-1) - ux*front_offset + vx*side_offset, ...
center_Y_list(end-1) - uy*front_offset + vy*side_offset];
B_point = [center_X_list(end-1) - ux*front_offset - vx*side_offset, ...
center_Y_list(end-1) - uy*front_offset - vy*side_offset];
C_point = [L_X + ux*front_offset - vx*side_offset, ...
L_Y + uy*front_offset - vy*side_offset];
D_point = [L_X + ux*front_offset + vx*side_offset, ...
L_Y + uy*front_offset + vy*side_offset];
rectangle_X = [A_point(1), B_point(1), C_point(1), D_point(1), A_point(1)];
rectangle_Y = [A_point(2), B_point(2), C_point(2), D_point(2), A_point(2)];
plot(rectangle_X, rectangle_Y, 'b-', 'LineWidth', 1.2);
if i > 1
bodyRects{end+1} = polyshape(rectangle_X, rectangle_Y);
else
head_front_pts = [A_point; B_point];
end
X = solution;
end
plot(center_X_list, center_Y_list, 'r-', 'LineWidth', 1.5);
for k = 1:length(bodyRects)
if any(isinterior(bodyRects{k}, head_front_pts(:,1), head_front_pts(:,2)))
disp(['发生碰撞,时间 t = ', num2str(t, '%.6f'), ' 秒']);
plot(head_front_pts(:,1), head_front_pts(:,2), 'ro', 'MarkerFaceColor','r');
collision_detected = true;
break;
end
end
drawnow;
if collision_detected
break;
end
end
第三问
clc;
clear;
close all;
A = (0.4503/(2*pi))^2;
R = @(theta) 0.4503/(2*pi)*theta;
numBody = 22;
head_theta0 = 32*pi;
syms theta
T = theta*sqrt(theta^2+1) + log(sqrt(theta^2+1) + theta);
C = subs(T, theta, head_theta0);
front_offset = 0.275;
side_offset = 0.15;
theta_vals = linspace(0, head_theta0, 2000);
spiral_X = R(theta_vals).*cos(theta_vals);
spiral_Y = R(theta_vals).*sin(theta_vals);
figure;
axis equal;
grid on;
xlabel('X'); ylabel('Y');
title('龙身体随时间运动的动画');
collision_detected = false;
for t = 0:10:400
clf; hold on; axis equal; grid on;
xlabel('X'); ylabel('Y');
title(['时间 t = ', num2str(t, '%.7f'), ' 秒']);
plot(spiral_X, spiral_Y, 'k--', 'LineWidth', 1);
question = T == C - 4*pi/0.4503*t;
head_theta = double(vpasolve(question, theta, [0, head_theta0]));
head_X = double(R(head_theta) * cos(head_theta));
head_Y = double(R(head_theta) * sin(head_theta));
center_X_list = head_X;
center_Y_list = head_Y;
X = head_theta;
plot(head_X, head_Y, 'ro', 'MarkerFaceColor','r');
ans = sqrt(head_Y ^ 2 + head_X ^ 2);
bodyRects = {};
for i = 1:numBody
if i == 1
L = 2.86; rect_length = 3.41;
else
L = 1.65; rect_length = 2.2;
end
syms Y
eqn = A*(Y^2 + X^2) - L^2 == 2*A*Y*X*cos(Y - X);
solution = double(vpasolve(eqn, Y, [X, X+pi/2]));
L_X = R(solution) * cos(solution);
L_Y = R(solution) * sin(solution);
plot(L_X, L_Y, 'bo', 'MarkerFaceColor','b');
center_X_list = [center_X_list, L_X];
center_Y_list = [center_Y_list, L_Y];
dx = L_X - center_X_list(end-1);
dy = L_Y - center_Y_list(end-1);
len = sqrt(dx^2 + dy^2);
ux = dx / len; uy = dy / len;
vx = -uy; vy = ux;
A_point = [center_X_list(end-1) - ux*front_offset + vx*side_offset, ...
center_Y_list(end-1) - uy*front_offset + vy*side_offset];
B_point = [center_X_list(end-1) - ux*front_offset - vx*side_offset, ...
center_Y_list(end-1) - uy*front_offset - vy*side_offset];
C_point = [L_X + ux*front_offset - vx*side_offset, ...
L_Y + uy*front_offset - vy*side_offset];
D_point = [L_X + ux*front_offset + vx*side_offset, ...
L_Y + uy*front_offset + vy*side_offset];
rectangle_X = [A_point(1), B_point(1), C_point(1), D_point(1), A_point(1)];
rectangle_Y = [A_point(2), B_point(2), C_point(2), D_point(2), A_point(2)];
plot(rectangle_X, rectangle_Y, 'b-', 'LineWidth', 1.2);
if i > 1
bodyRects{end+1} = polyshape(rectangle_X, rectangle_Y);
else
head_front_pts = [A_point; B_point];
end
X = solution;
end
plot(center_X_list, center_Y_list, 'r-', 'LineWidth', 1.5);
for k = 1:length(bodyRects)
if any(isinterior(bodyRects{k}, head_front_pts(:,1), head_front_pts(:,2)))
disp(['发生碰撞,时间 t = ', num2str(t, '%.6f'), ' 秒']);
plot(head_front_pts(:,1), head_front_pts(:,2), 'ro', 'MarkerFaceColor','r');
collision_detected = true;
break;
end
end
drawnow;
if collision_detected
break;
end
end
fprintf('%.15f',ans);