[原创]用clough书中的逐步积分法求非线性系统数值解的源程序

scunwzh

利用逐步积分法求duffning系统的解

如果不怕修改参数麻烦的话采用采用如下程序（matlab6.5），此种方法运行速度快。

% y''+2y'+10^4y+5*10^9y^3=4*cos(30*t^2)
% define the preload
clc;
clear all;
close all;
deltat=0.002;
t=0:deltat:6;
p=4*cos(30*t.^2);
% initial the displacement and velocity
y=0;v=1;m=1;c=2;
len=length(t)-1;
displacement=zeros([1,len]);
velocity=zeros([1,len]);
accelaration=zeros([1,len]);
for i=1:len
k=10^4+15*10^9*y^2;
a=(p(i)-c*v-10^4*y-5*10^9*y^3)/m;
equivelentk=k+6*m/(deltat^2)+3/deltat*c;
deltap=(p(i+1)-p(i))+m*(6/deltat*v+3*a)+c*(3*v+deltat/2*a);
deltay=deltap/equivelentk;
deltav=3/deltat*deltay-3*v-deltat/2*a;
y=y+deltay;
v=v+deltav;
displacement(i)=y;
velocity(i)=v;
accelaration(i)=a;
end
figure(1)
plot(t(1:len),displacement)
figure(2)
plot(t(1:len),velocity)
figure(3)
plot(t(1:len),accelaration)

如果懒得每次修改参数可采用如下程序，程序中使用了symbolic object，运行速度慢

function [displacement,velocity,accelaration]=intdyn(fc,fk,pt,m,x0,v0,deltat,tt)
% fc ,fk,pt are symbolic object,fc is damping force,fk is elastic force,pt
% is preload.
% fc=fc(v);fk=fk(x);pt=pt(t);
% pt is numeric array of preload
% deltat is length of time step
% x0 and v0 are the initial data of displacement and velocity,respectly
% tt is time of sampling
t1=0:deltat:tt;
len=length(t1)-1;
for i=1:len+1
p(i)=subs(pt,'t',t1(i));
end
% displacement=zeros([1,len]);
% velocity=zeros([1,len]);
% accelaration=zeros([1,len]);
xx=x0;
vv=v0;
k1=diff(fk,'x');
c1=diff(fc,'v');
for i=1:len
% replace the symbolic object variable
k=subs(k1,'x',xx);
c=subs(c1,'v',vv);
ffk=subs(fk,'x',xx);
ffc=subs(fc,'v',vv);
% decrease the accumulate of error
aa=(p(i)-ffc-ffk)/m;
% compute the equivelent stiff
equivelentk=k+6*m/(deltat^2)+3/deltat*c;
% compute the equivelent load
deltap=(p(i+1)-p(i))+m*(6/deltat*vv+3*aa)+c*(3*vv+deltat/2*aa);
% compute the increment of displacement
deltax=deltap/equivelentk;
% compute the increment of velocity
deltav=3/deltat*deltax-3*vv-deltat/2*aa;
% compute the displacement and velocity,respectly
xx=xx+deltax;
vv=vv+deltav;
% output
displacement(i)=xx;
velocity(i)=vv;
accelaration(i)=aa;
end
采用如下程序调用以上子程序

clc;clear all;close all
syms x v t;
deltat=0.002;
tt=6;
x0=0;v0=1;m=1;
[dx,dv,da]=intdyn(2*v,10^4*x+5*10^9*x^3,4*cos(30*t.^2),m,x0,v0,deltat,tt);

初次发帖。不对之处请指正。

gghhjj

好东西，非常感谢，明天试试看

linqus

支持原创！！

希望大家对程序进行测试以资楼主改进。、
之后设置为精华贴。

[此贴子已经被作者于2005-11-9 23:30:04编辑过]

scunwzh

上面文件是clough书中的内容，有助于理解上述程序。我总结的与clough的不同之处在于将瞬时位移和瞬时速度展开为Taylor级数，然后利用线性加速度的概念回代。这样比clough书和其它的一些数值计算的书要好理解一些。