maker系列－强迫振动欧拉数值解法程序

MVH

maker系列－强迫振动欧拉数值解法程序
function maker3(m,c,k,x0,v0,time,w,f0,delt)
%oule改进法，利用泰勒展开式
%m-为质量，c-阻尼，k-刚度，x0-初位移，v0-除速度；time-仿真时间，delt离散时间
%mx(t)''+cx(t)'+kx(t)=f(t)
%x(t)''=(f(t)-kx(t)-cx(t)'))/m
%x(t)'''=(f(t)'-kx(t)'-cx(t)'')/m
%x(t+delt)'=x(t)'+x(t)''*delt+x(t)'''*delt*delt/2，可以用Simpson或梯形公式求得，会更加精确
%x(t+delt)=x(t)+x(t)'*delt+x(t)''*delt*delt/2
%chinamaker@dytrol.com, 动力学与控制论谈:http://www.dytrol.com
%2003.11.27

i=1;
n=time/delt;
disp=zeros(n,1);
for t=0:delt:time
xdd=(f0*sin(w*t)-k*x0-c*v0)/m;
xddd=(f0*w*cos(w*t)-k*v0-c*xdd)/m;
xd=v0+xdd*delt+xddd*delt*delt/2;
x=x0+xd*delt+xdd*delt*delt/2;
disp(i)=x;
i=i+1;
x0=x;
v0=xd;
end
t=0:delt:time;
plot(t,disp);
grid;
xlabel('time(s)');
ylabel('Amplitude');

实例计算结果:
>> maker3(20,1.5,40,1.5,1.5,100,5,100,0.0001)

MVH

function maker4(m,c,k,x0,v0,time,delt)
%线性加速度法计算三自由度简谐受迫振动:改编于尚涛等的书
%m-为质量，c-阻尼，k-刚度，x0-初位移，v0-除速度；time-仿真时间，delt-时间步长
% MX''+CX'+KX=F
% X(t+delt)=[M+delt*C/2+delt^2*K/6]^(-1){F(t+delt)-C[X(t)+delt*X(t)'']-K[X(t)+delt*X(t)'+delt^2*X(t)''/3]}
% X(t+delt)'=X(t)'+delt*[X(t)''+X(t+delt)]/2
% X(t+delt)=X(t)+delt*X(t)'+(delt^2)*X(t)''/3+(delt^2)*X(t+delt)''/6
%chinamaker@dytrol.com
%2003.11.27

n=time/delt;
disp=zeros(3,n);%设定存储位移矩阵
m_inv=inv(m+delt*c/2+delt^2*k/6);
[mod,fre]=eig(inv(m)*k);
diag(sqrt(fre));%固有频率
i=1;
for t=0:delt:time
f=[2.0*sin(3.5*t) -2.0*cos(2*t) 1.0*sin(3*t)]';%外扰力
if t==0
xdd0=inv(m)*(f-k*x0-c*v0);%初始加速度
else
xdd=m_inv*(f-c*(v0+delt/2*xdd0)-k*(x0+delt*v0+(1/2-1/6)*delt^2*xdd0));%计算加速度
x=m_inv*(m*(x0+delt*v0+delt^2/3*xdd0)+c*(delt/2*x0+delt^2/3*v0+delt^3/12*xdd0)+delt^2/6*f);%计算位移
xd=v0+delt*(xdd0+xdd)/2;%计算速度
xdd0=xdd;
v0=xd;
x0=x;
disp(:,i)=x0;
i=i+1;
end
end
t=1:n;
figure('numbertitle','off','name','自由度1的位移','pos',[450 180 400 420]);
plot(t,disp(1,:)),grid,xlabel('时间(s/10)'),title('自由度1的时程曲线');
figure('numbertitle','off','name','自由度2的位移','pos',[450 180 400 420]);
plot(t,disp(2,:)),grid,xlabel('时间(s/10)'),title('自由度2的时程曲线');
figure('numbertitle','off','name','自由度3的位移','pos',[450 180 400 420]);
plot(t,disp(3,:)),grid,xlabel('时间(s/10)'),title('自由度3的时程曲线');
%end

实例：
m=2*[1 0 0;0 1 0;0 0 1];
c=1.5*[2 -1 0;-1 2 -1;0 -1 2];
k=50*[2 -1 0;-1 2 -2;0 -1 2];
x0=[1 1 1]';
v0=[1 1 1]';
delt=0.1;
time=100;
作用力由于用到较多的数据，直接写在原程序中。

>> maker4(m,c,k,x0,v0,time,delt)
结果如下

MVH

maker系列-超松弛迭代法
考试练手写的简单程序,不要见笑,呵呵.

function x=sor(A,b,w)
%超松弛迭代法
%A*x=b
%w为松弛系数，1<=w<=2;
numer=length(A);
D=zeros(numer);
L=zeros(numer);
U=zeros(numer);
x0=zeros(numer,1);
for i=1:numer
D(i,i)=A(i,i);%对角矩阵
x0(i,1)=1;%初始迭代向量
end
for i=2:3
for j=1:i-1
U(j,i)=-A(j,i);%U上半部矩阵
end
for j=1:i-1
L(i,j)=-A(i,j); %L下半部矩阵
end
end
Lw=inv(D-w*L)*((1-w)*D+w*U);%迭代矩阵
g=inv(D-w*L)*w*b;
errs=1;
while errs>0.001
x1=Lw*x0+g;
errs=norm(x1-x0);
x0=x1;
end
x=x0;

实例
>> A=[4,3,0;3,4,1;0,-1,4];
>> b=[24;20;-24];
>> w=1.25;
>> x=sor(A,b,w)

x =

3.0001
4.0000
-5.0000

dingxinran

MVH 发表于 2005-7-28 17:47

                               
登录/注册后可看大图

function maker4(m,c,k,x0,v0,time,delt)
%线性加速度法计算三自由度简谐受迫振动:改编于尚涛等的书
%m-为 ...

>> m=2*[1 0 0;0 1 0;0 0 1];
>> c=1.5*[2 -1 0;-1 2 -1;0 -1 2];
>> k=50*[2 -1 0;-1 2 -2;0 -1 2];
>> x0=[1 1 1]';
>> v0=[1 1 1]';
>> delt=0.1;
>> time=100;
>> maker4(m,c,k,x0,v0,time,delt)
??? Input argument "delt" is undefined.

Error in ==> maker4 at 9
n=time/delt;

运行之后出现这个错误是什么意思呢 怎么改呢