声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

声振论坛 展示 振动理论 非线性振动 查看内容

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

2016-3-21 09:50| 发布者: aspen| 查看: 1619| 评论: 3|原作者: scunwzh|来自: 声振论坛

摘要: 利用逐步积分法求duffning系统的解如果不怕修改参数麻烦的话采用采用如下程序(matlab6.5),此种方法运行速度快。
利用逐步积分法求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 2005-11-9 22:34
好东西,非常感谢,明天试试看
引用 linqus 2005-11-9 23:25
支持原创!!

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

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

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

查看全部评论(3)

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-5-10 04:53 , Processed in 0.050238 second(s), 23 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

返回顶部