MATLAB程序设计课后作业五

一、改写高斯积分函数

clc;clear;close all;
ff=@(x) x*x*exp(x); 
s1 = GaussL(ff,0,1,2);% 积分点为2点
s2 = GaussL(ff,0.5,1,3); % 积分点为3点
disp(s1);disp(s2);

function s = GaussL(f,a,b,n)
% 输入：f为待积分函数的函数句柄，a和b为积分上下限，n为积分点个数
% 输出：s为积分值
if n == 2
    x = [-0.577350269189626, 0.577350269189626]; % 高斯点
    w = [1, 1]; % 权重
elseif n == 3
    x = [-0.774596669241483, 0, 0.774596669241483]; % 高斯点
    w = [0.555555555555556, 0.888888888888889, 0.555555555555556]; % 权重
else
    error("n必须为2或者3"); 
end
% 计算积分值
s = 0;
for i = 1:n
    s = s + w(i) * feval(f, (b-a)/2*x(i) + (b+a)/2); 
end
s = s * (b-a)/2; 
end

二、处理污水

clc;clear;close all;
% 定义优化变量
x = optimvar('x', 'LowerBound', 0, 'UpperBound', 2e4);
y = optimvar('y', 'LowerBound', 0, 'UpperBound', 1.4e4);
prob = optimproblem('ObjectiveSense', 'minimize');
% 定义目标函数
prob.Objective = 0.1*x + 0.08*y;
% 添加约束条件
prob.Constraints.cons = (2e4 - x)*0.8 + 1.4e4 - y <= 0.002*(5e6 + 2e4 - x + 2e5);
% 求解优化问题
sol = solve(prob);
% 显示最优解和目标函数值
disp('第一个化工厂每天应处理的污水量(立方米):')
disp(sol.x);
disp('第二个化工厂每天应处理的污水量(立方米):')
disp(sol.y);
disp('两厂总的处理污水费用(元):')
disp(prob.Objective.evaluate(sol));

三、各种梁受均布载荷的静定问题

clc;clear;close all;
% 定义求解参数和区间
L = 10; 
q = 1000;
E = 2e11; 
I = 1e-4; 
x = linspace(0, L, 100); % 区间
solinit = bvpinit(x, [0 0 0 0]); % 初值
% 调用bvp4c函数求解，并画出挠度、转角和弯矩图
type = '悬臂梁';
sol = bvp4c(@(x, y) odefun(x, y, q, E, I), @(ya, yb) bcfun(ya, yb, type), solinit); % 求解
y = sol.y; 
w = y(1, :); 
theta = y(2, :);
M = -E*I*y(3, :);
subplot(3, 1, 1)
plot(x, w, 'b')
xlabel('x')
ylabel('w')
title(['挠度图（' type '）'])
subplot(3, 1, 2)
plot(x, theta, 'r')
xlabel('x')
ylabel('\theta')
title(['转角图（' type '）'])
subplot(3, 1, 3)
plot(x, M, 'g')
xlabel('x')
ylabel('M')
title(['弯矩图（' type '）'])

type后边的参数需要自行修改，分别填入两端固支梁，两端简支梁和悬臂梁后分别运行得结果！！！

function bc = bcfun(ya, yb, type)
switch type
    case '两端简支梁'
        bc = [ya(1); yb(1); ya(3); yb(3)];
    case '两端固支梁'
        bc = [ya(1); yb(1); ya(2); yb(2)];
    case '悬臂梁'
        bc = [ya(1); ya(2); ya(3); ya(4)];
end
end

% 定义梁的边界条件和微分方程
function dydx = odefun(x, y, q, E, I)
dydx = [y(2); y(3); y(4); -q/(E*I)];
end

四、外伸梁的拓展

clc;clear;close all;
% 定义外神梁的参数
L = 10; 
E = 2e11; 
I = 0.01; 
q = 1000; 


f = @(x,y) [y(2); y(3); y(4); q/(E*I)];
bc = @(ya,yb) [ya(1); ya(2); yb(1); yb(2)];

xspan = [0 L];
yinit = [0; 0; 0; 0];
solinit = bvpinit(xspan,yinit);

% 调用bvp4c函数求解边值问题
sol = bvp4c(f,bc,solinit);

% 提取解向量
x = sol.x; 
v = sol.y(1,:); 
theta = sol.y(2,:); 
M = -E*I*sol.y(3,:); 

% 绘制挠度，转角和弯矩图
figure(1)
plot(x,v,'b')
xlabel('x (m)')
ylabel('v (m)')
title('挠度图')
grid on

figure(2)
plot(x,theta,'r')
xlabel('x (m)')
ylabel('\theta (rad)')
title('转角图')
grid on

figure(3)
plot(x,M,'g')
xlabel('x (m)')
ylabel('M (N*m)')
title('弯矩图')
grid on

第五题见MATLAB程序设计课后作业六，别问我为什么分开发，问就是为了获得创作推广。