『求助』Newmark法的Matlab程序
浏览:128864
function [ua,ua1,ua2]=newmark(m,k,c,pt,u,u1,u2,alpha,dertat,derta,n)
%by swell2005
%--------
format long
a0=1/(alpha*dertat^2)
a1=derta/(alpha*dertat)
a2=1/(alpha*dertat)
a3=0.5/alpha-1
a4=derta/alpha-1
a5=dertat*(derta/alpha-2)/2
a6=dertat*(1-derta)
a7=derta*dertat
k=k+a0*m+a1*c
[r,s]=ldlt(k);
%-----------
for i=0:n
pta=pt+m*(a0*u+a2*u1+a3*u2)+c*(a1*u+a4*u1+a5*u2);
uaa=r\pta;
uaa=s\uaa;
ua(:,i+1)=r'\uaa;
ua2(:,i+1)=a0*(ua(:,i+1)-u)-a2*u1-a3*u2;
ua1(:,i+1)=u1+a6*u2+a7*ua2(:,i+1);
u=ua(:,i+1);u1=ua1(:,i+1);u2=ua2(:,i+1);
end
format
我编写的过程基本上跟以上程序段是一致的,只是其中的pt这个变量,我要求解的动力学方程里,pt是一个关于u的函数,而u本身有应该是一个关于时间t的函数,我在不知道u(t)的确切表达式的情况下,该如何编制这段程序?
%by swell2005
%--------
format long
a0=1/(alpha*dertat^2)
a1=derta/(alpha*dertat)
a2=1/(alpha*dertat)
a3=0.5/alpha-1
a4=derta/alpha-1
a5=dertat*(derta/alpha-2)/2
a6=dertat*(1-derta)
a7=derta*dertat
k=k+a0*m+a1*c
[r,s]=ldlt(k);
%-----------
for i=0:n
pta=pt+m*(a0*u+a2*u1+a3*u2)+c*(a1*u+a4*u1+a5*u2);
uaa=r\pta;
uaa=s\uaa;
ua(:,i+1)=r'\uaa;
ua2(:,i+1)=a0*(ua(:,i+1)-u)-a2*u1-a3*u2;
ua1(:,i+1)=u1+a6*u2+a7*ua2(:,i+1);
u=ua(:,i+1);u1=ua1(:,i+1);u2=ua2(:,i+1);
end
format
我编写的过程基本上跟以上程序段是一致的,只是其中的pt这个变量,我要求解的动力学方程里,pt是一个关于u的函数,而u本身有应该是一个关于时间t的函数,我在不知道u(t)的确切表达式的情况下,该如何编制这段程序?