姜启源老师数学模型课本相关主题MATALB代码，请自行验证一下

姜启源老师的数学模型课本，每个部分都有具体的实验内容和示例代码，以便学生更好地理解和应用这些数学模型。

1. 数学建模概述

实验内容：

• 目标：熟悉建模流程，通过“椅子问题”体验建模的基本思路。
• 任务：
• 描述椅子放置的最佳位置问题。
• 提出简化假设。
• 建立数学模型并求解。
• 结果解释与讨论。

% 椅子问题示例：在一个矩形房间内找到最佳放置两张椅子的位置
room_width = 5; % 房间宽度
room_length = 4; % 房间长度
chair_radius = 0.5; % 椅子半径

% 定义目标函数：最大化两把椅子中心的距离
objective_function = @(positions) -sqrt((positions(1)-positions(3))^2 + (positions(2)-positions(4))^2);

% 设置边界条件
lb = repmat(chair_radius, 1, 4); % 下界
ub = repmat(room_width-chair_radius, 1, 4); % 上界

initial_guess = [1, 1, room_width-chair_radius, room_length-chair_radius]; % 初始猜测

options = optimoptions(@fmincon, 'Display', 'iter');

[best_positions, min_distance] = fmincon(objective_function, initial_guess, [], [], [], [], lb, ub, [], options);

disp(['Best positions for chairs:', best_positions']);
disp(['Minimum distance between chair centers: ', num2str(-min_distance)]);

2. 初等模型

实验内容：

• 目标：建立比例、代数、几何模型的思维框架，通过“资源分配问题”体验简单建模的实际应用。
• 任务：
• 定义资源分配问题。
• 建立方程组并求解。
• 结果解释与讨论。

% 资源分配问题示例：给定预算分配两种商品的数量以最大化利润
budget = 1000; % 总预算
cost_productA = 100; % 商品A的成本单价
profit_productA = 50; % 商品A的利润单价
cost_productB = 200; % 商品B的成本单价
profit_productB = 80; % 商品B的利润单价

% 定义目标函数：最大化总利润
objective_function = @(quantities) -(quantities(1)*profit_productA + quantities(2)*profit_productB);

% 设置约束条件
constraints_matrix = [
    cost_productA, cost_productB; % 成本不超过预算
];

rhs_constraints = budget;

lower_bounds = zeros(2, 1); % 数量不低于0

options = optimoptions('@linprog', 'Display', 'final');

[opt_quantities, opt_profit] = linprog(fvec, Aineq, bineq, Aeq, beq, lower_bounds, upper_bounds, options);

opt_profit = -opt_profit; % 将负号去掉恢复正值

disp(['Optimal quantity of Product A: ', num2str(opt_quantities(1))]);
disp(['Optimal quantity of Product B: ', num2str(opt_quantities(2))]);
disp(['Maximum profit: ', num2str(opt_profit)]);

3. 优化模型

实验内容：

• 目标：认识线性规划与单纯形法的基本原理，通过“生产计划问题”建立优化思维并实践编程（MATLAB/Python）。
• 任务：
• 定义生产计划问题。
• 使用 linprog 函数求解线性规划问题。
• 结果解释与讨论。

% 生产计划问题示例：最大化产量同时满足原材料限制
maximize_production = [-1, -2]; % 目标系数（取反因为默认是最小化）

material_limits = [
    2, 1; % 材料A限制
    1, 2; % 材料B限制
];

available_materials = [100, 100]; % 可用材料总量

product_requirements_lower_bound = [0, 0]; % 产品数量最低为0

options = optimoptions('linprog', 'Algorithm', 'simplex-largescale', 'Display', 'none');

[solution, optimal_value, exitflag] = linprog(maximize_production, material_limits, available_materials, [], [], product_requirements_lower_bound, []);

if exitflag > 0
    disp(['Optimal production plan:'])
    disp(solution)
    disp(['Total production value maximized: ', num2str(-optimal_value)])
else
    disp('The problem does not have a feasible solution.')
end

4. 微分方程模型

实验内容：

• 目标：认识常微分方程建模的基本方法，通过“SIR传染病模型”建立动态系统建模思维并完成求解与分析。
• 任务：
• 构建 SIR 模型。
• 使用 ode45 函数求解 ODE 方程组。
• 结果可视化与讨论。

function sirModel()
    % 参数设置
    beta = 0.3; % 传染率
    gamma = 0.1; % 康复率
    
    % 初始化人口状态
    N = 1000; % 总人数
    I0 = 1; % 初始感染人数
    R0 = 0; % 初始康复人数
    S0 = N - I0 - R0; % 初始易感人数
    
    tspan = [0 150]; % 时间区间
    
    % 初始条件向量
    y0 = [S0, I0, R0];
    
    % 求解ODE
    [~, Y] = ode45(@(t,y) sirdynamics(beta, gamma, y), tspan, y0);
    
    % 数据提取
    S = Y(:,1);
    I = Y(:,2);
    R = Y(:,3);
    
    % 绘制结果
    figure;
    plot(tspan, S/N, '-b', 'LineWidth', 2);
    hold on;
    plot(tspan, I/N, '-r', 'LineWidth', 2);
    plot(tspan, R/N, '-g', 'LineWidth', 2);
    legend('Susceptible', 'Infectious', 'Recovered');
    xlabel('Days');
    ylabel('Fraction of Population');
    title('SIR Epidemic Model Dynamics');
    grid on;
end

function dydt = sirdynamics(beta, gamma, y)
    S = y(1);
    I = y(2);
    R = y(3);
    N = S + I + R;
    
    dSdt = -beta*S*I/N;
    dIdt = beta*S*I/N - gamma*I;
    dRdt = gamma*I;
    
    dydt = [dSdt; dIdt; dRdt];
end

sirModel();

5. 差分方程模型

实验内容：

• 目标：认识离散模型与稳定性分析的基本原理，通过“种群动态模型”建立离散系统建模思维。
• 任务：
• 构建 Logistic 种群增长模型。
• 使用循环迭代求解差分方程。
• 结果可视化与讨论。

% Logisitic Growth Model Example
clear all;
close all;

% Parameters
alpha = 0.2; % 自然增长率
rho = 0.01; % 密度依赖因子
N_max = 1000; % 最大人口容量
time_steps = 100; % 时间步长
population_0 = 100; % 初始人口

% Initialize population array
populations = zeros(time_steps, 1);
populations(1) = population_0;

% Iterate through each time step
for t = 1:(time_steps-1)
    populations(t+1) = populations(t) + alpha * populations(t) * (1 - rho * populations(t) / N_max);
end

% Plot the result
figure;
plot(populations, 'LineWidth', 2);
xlabel('Time Steps');
ylabel('Population Size');
title('Logistic Growth Model');
grid on;

6. 概率模型

实验内容：

• 目标：认识随机过程与马尔可夫链的基本概念，通过“天气预测模型”建立概率建模思维并完成设计与验证。
• 任务：
• 构建 Markov Chain 天气预测模型。
• 使用转移矩阵进行长期预测。
• 结果可视化与讨论。

% Weather Prediction Model using Markov Chains
states = {'Sunshine', 'Rain'};
transition_probabilities = [
    0.9, 0.1; % P(Sun -> Sun), P(Sun -> Rain)
    0.5, 0.5  % P(Rain -> Sun), P(Rain -> Rain)
];

current_state = states{1}; % Initial state: Sunshine
days_to_predict = 100;

weather_sequence = cell(days_to_predict, 1);
weather

7. 统计模型MATLAB实验建议

实验内容：

• 目标：理解回归分析与数据拟合的基本方法，通过“房价预测”建立数据驱动建模思维并实践（Python/R）。
• 任务：
• 使用线性回归模型分析房屋特征与价格之间的关系。
• 收集或使用示例房地产数据集。
• 在MATLAB中实现线性回归或使用Statistics and Machine Learning Toolbox。

% 假设已有数据集houses.csv，包含房屋面积sqft和价格price
data = readtable('houses.csv');
X = data(:, 'sqft'); % 特征
Y = data(:, 'price'); % 目标变量

% 分割数据为训练集和测试集
cv = cvpartition(height(X), 'HoldOut', 0.3);
Xtrain = X(training(cv), :);
Ytrain = Y(training(cv), :);
Xtest = X(test(cv), :);
Ytest = Y(test(cv), :);

% 线性回归模型
mdl = fitlm(Xtrain, Ytrain);

% 预测
Ypred = predict(mdl, Xtest);

% 评估模型
disp(['R-squared: ', num2str(mdl.Rsquared.Ordinary)]);

8. 图论模型MATLAB实验建议

实验内容：

• 目标：认识最短路径与最小生成树算法的基本原理，通过“交通网络优化”建立图论建模思维并实践编程（Python）。
• 任务：
• 构建一个交通网络图，定义节点和边的权重。
• 实现Dijkstra算法或使用MATLAB内置函数找到两节点间的最短路径。
• 应用Prim's或Kruskal's算法找到最小生成树。

% 假设交通网络数据
G = graph(...
    [1 2 2 3 3 4 4 5], ... % 边的起点
    [2 3 4 4 5 5 6 6], ... % 边的终点
    [10 5 15 12 8 7 9 11]  % 边的权重
);

% 最短路径
sourceNode = 1;
targetNode = 6;
[shortestPath, dist] = shortestpath(G, sourceNode, targetNode);
fprintf('Shortest path from %d to %d: %d, Distance: %d\n', sourceNode, targetNode, shortestPath, dist);

% 最小生成树
if strcmp(version, '9.10') % 假定MATLAB版本支持最小生成树函数
    mst = minspantree(G);
    edges = edges(mst);
    fprintf('Edges in the minimum spanning tree: %d\n', edges);
else
    disp('Upgrade MATLAB for built-in minspantree function or implement manually.');
end

9. 非线性规划

实验内容：

• 目标：学会使用 fmincon 函数解决非线性规划问题。
• 任务：
• 定义一个简单的非线性目标函数和约束条件。
• 使用 fmincon 找到最优解。
• 绘制目标函数曲面和约束区域，可视化结果。

% 示例：最小化 f(x,y) = (x-3)^2 + y^2 主观约束 g(x,y) ≤ 0
fun = @(x) (x(1)-3).^2 + x(2).^2;
nonlcon = @(x)[-(x(1)+x(2)); -(x(1)*x(2))];
options = optimoptions('fmincon','Display','iter');
[x,fval] = fmincon(fun,[0;0],[],[],[],[],[-Inf,-Inf],[Inf,Inf],nonlcon,options);

% 可视化
[X,Y] = meshgrid(-5:0.1:5,-5:0.1:5);
Z = fun([X(:), Y(:)]).';
contour(X,Y,Z,logspace(-2,2,20));
hold on;
plot([-Inf Inf],[-Inf Inf],'r--'); % 约束区域
scatter(x(1),x(2),'filled')
title('Nonlinear Optimization with Constraints')
xlabel('x')
ylabel('y')
legend('Contour of Objective Function','Constraints','Optimal Point')

10. 动态规划

实验内容：

• 目标：理解动态规划的基本原理并通过背包问题实例演示。
• 任务：
• 编写递归和迭代版本的动态规划算法。
• 解决经典的0-1背包问题。
• 计算时间和空间复杂度。

% 示例：0-1 Knapsack Problem using Dynamic Programming
W = 50; % Maximum weight capacity
weights = [10, 20, 30];
values = [60, 100, 120];
n = length(weights);

dp = zeros(n+1,W+1);

for i=1:n
    for w=1:W
        if weights(i) > w
            dp(i+1,w) = dp(i,w);
        else
            dp(i+1,w) = max(dp(i,w), values(i) + dp(i,max(w-weights(i),0)));
        end
    end
end

disp(['Maximum value in knapsack is ', num2str(dp(end,end))]);

11. 排队论模型

实验内容：

• 目标：模拟M/M/1排队系统。
• 任务：
• 设计参数（到达率λ和服务率μ）。
• 进行离散事件模拟。
• 分析平均等待时间、队长分布等指标。

% M/M/1 Queue Simulation Example
lambda = 1; % Arrival rate
mu = 2;     % Service rate
T_sim = 1000; % Total simulation time

arrival_times = exprnd(1/lambda,T_sim,1);
service_times = exprnd(1/mu,sum(arrival_times < T_sim));

queue_time = cumsum(service_times) - arrival_times(find(arrival_times < T_sim));
waiting_time = queue_time - service_times;

mean_waiting_time = mean(waiting_time(queue_time >= 0))
histogram(waiting_time(queue_time >= 0), 'Normalization', 'pdf')

figure;
subplot(2,1,1);
stairs(sort(queue_time(queue_time>=0)), linspace(0, sum(queue_time(queue_time>=0)>0)/length(queue_time(queue_time>=0)), length(queue_time(queue_time>=0)))
title('Queue Time Distribution')
xlabel('Time')
ylabel('Probability Density')

subplot(2,1,2)
stairs(sort(waiting_time(queue_time>=0)), linspace(0, sum(waiting_time(queue_time>=0)>0)/length(waiting_time(queue_time>=0)), length(waiting_time(queue_time>=0)))
title('Waiting Time Distribution')
xlabel('Time')
ylabel('Probability Density')

12. 博弈论模型

实验内容：

• 目标：计算纳什均衡。
• 任务：
• 构造支付矩阵。
• 使用 game-theory-toolbox 或自定义函数寻找纯策略纳什均衡。

% Prisoner's Dilemma Game Matrix
payoff_A = [-1 -3; 0 -2];
payoff_B = payoff_A'; % Symmetric game

% Find Nash Equilibrium
NE_toolbox = gametoolbox(payoff_A, payoff_B);
ne_pure_strategies = NE_toolbox.puresolve();

if ~isempty(ne_pure_strategies)
    disp('Pure Strategy Nash Equilibria:');
    ne_pure_strategies'
else
    disp('No pure strategy Nash equilibria found.');
end

注意：gametoolbox 是第三方工具箱，可能需要安装。

13. 网络流模型

实验内容：

• 目标：求解最大流问题。
• 任务：
• 创建网络图。
• 使用 mldivide 或 graph 工具箱求解最大流。

s = 1;      % Source node
t = 5;      % Sink node
capacity = sparse([
    1 2 16;
    1 3 13;
    2 3 10;
    2 4 12;
    3 2 4;
    3 4 9;
    3 5 14;
    4 3 6;
    4 5 7]);

G = digraph(sparse(capacity(:,1), capacity(:,2), ones(size(capacity,1),1)));

flowMatrix = flow(G,s,t,capacity(:,3));

total_flow = sum(flowMatrix(s,:))

% Visualize the network and flows
h = plot(G,'EdgeLabel',num2cell(string(flowMatrix.'./full(capacity))));
highlight(h,G.Lineage(t),"Color","red")
title("Max Flow from Node ",string(s),"to Node ", string(t))

13. 随机模拟

实验内容：

• 目标：使用蒙特卡洛方法估计π值。
• 任务：
• 投掷大量的随机点到单位圆内。
• 计算落在圆内的比例。
• 估算 π 的近似值。

% Monte Carlo Method to Estimate Pi
N = 1e6; % Number of random points
points_x = rand(N,1);
points_y = rand(N,1);

inside_circle = sqrt(points_x.^2 + points_y.^2) <= 1;
pi_estimate = 4 * sum(inside_circle) / N;

fprintf('Estimated pi: %.6f\n', pi_estimate)

% Visualization
in_points = inside_circle == true;
out_points = inside_circle == false;

figure;
plot(points_x(in_points), points_y(in_points), '.', 'MarkerSize', 1, 'DisplayName', 'Inside Circle');
hold on;
plot(points_x(out_points), points_y(out_points), '.', 'MarkerSize', 1, 'DisplayName', 'Outside Circle');
axis equal;
xlim([0 1])
ylim([0 1])
title(sprintf('Monte Carlo Estimation of Pi (%d Points)', N))
legend show

14. 机器学习模型

实验内容：

• 目标：执行主成分分析（PCA）。
• 任务：
• 加载数据集。
• 应用 PCA 进行特征提取。
• 视觉化降维效果。

load fisheriris.mat
data = meas;
labels = species;

pca_model = fitzscore(data);
coefficients = pca_model.Coefficients;
scores = data*pca_model.Scores';

% Plotting first two principal components
gscatter(scores(:,1), scores(:,2), labels, [], '.')
title('First Two Principal Components')
xlabel('Principal Component 1')
ylabel('Principal Component 2')
legend unique(labels)