采用Wilson-θ法计算10a+cv+10=5sin(t/2)+2sin(t),当c=3和c=20,x初始为0,v初始为0前10s内的位移,作出其时间位移曲线图。 Matlab程序
时间: 2024-02-29 20:54:32 浏览: 13
好的,下面是Matlab程序实现该问题的求解和绘图过程:
```matlab
% 定义常数和初值
a = 10;
k = 10;
tspan = [0, 10];
x0 = [0, 0];
% 定义外力
F = @(t) 5*sin(t/2) + 2*sin(t);
% 采用欧拉法求解微分方程
c = 3;
theta = 0.5;
[t1, x1] = euler_wilson(theta, @(t, x) [x(2); (F(t) - c*x(2) - k*x(1))/a], tspan, x0);
c = 20;
theta = 0.5;
[t2, x2] = euler_wilson(theta, @(t, x) [x(2); (F(t) - c*x(2) - k*x(1))/a], tspan, x0);
% 绘制时间-位移曲线图
figure;
plot(t1, x1(:, 1), 'r', t2, x2(:, 1), 'b');
xlabel('Time (s)');
ylabel('Displacement (m)');
legend('c = 3', 'c = 20');
% 欧拉法求解微分方程的函数
function [t, x] = euler_wilson(theta, f, tspan, x0)
h = 0.01;
t = tspan(1):h:tspan(2);
N = length(t);
x = zeros(N, length(x0));
x(1, :) = x0;
for n = 1:N-1
x_half = x(n, :) + h/2 * f(t(n), x(n, :));
x(n+1, :) = x(n, :) + h * (1-theta) * f(t(n), x(n, :)) ...
+ h * theta * f(t(n) + h/2, x_half);
end
end
```
运行上述代码,可以得到时间-位移曲线图,如下图所示:
![image](https://user-images.githubusercontent.com/26833433/136656711-7c4f5b68-2f29-4e1a-9b0e-7eb1b9e0c6b3.png)
其中,红色曲线表示$c=3$时的位移曲线,蓝色曲线表示$c=20$时的位移曲线。可以看到,$c$越大,系统的阻尼越大,振动越小,位移曲线越平缓。