用的是wilson-θ方法,M=[1 0;0 1],K=[8 -2;-2 4],C=[0.27 -0.04;-0.04 0.19],P=[0;sint],θ=1.4,时间步长0.1s。
请问下这个p该怎么输入,程序如下
function [u u1 u2]=multi_wiltheta(M,K,C,d,e,f,P,theta,t)
n=length(P(:,1));m=length(P(1,:));
u=zeros(n,m);u1=zeros(n,m);u2=zeros(n,m);uu=zeros(n,m);PP=zeros(n,m);
u(:,1)=d;u1(:,1)=e;u2(:,1)=f;
l=6/(theta*t)^2;q=3/(theta*t);r=6/(theta*t);s=theta*t/2;
KK=K+l*M+q*C;
for i=1:m-1
pp(:,i)=P(:,i)+theta*(P(:,i+1)-P(:,i))+M*(l*u(:,i)+r*u1(:,i)+2*u2(:,i))+C*(q*u(:,i)+2*u1(:,i)+s*u2(:,i));
uu(:,i)=KK\pp(:,i);
u2(:,i+1)=u2(:,i)+l/theta*(uu(:,i)-u(:,i))-l*t*u1(:,i)-3/theta*u2(:,i);
u1(:,i+1)=u1(:,i)+t/2*(u2(:,i+1)+u2(:,i));
u(:,i+1)=u(:,i)+t*u1(:,i)+t^2/6*(u2(:,i+1)+2*u2(:,i));
end
end
请问下这个p该怎么输入,程序如下
function [u u1 u2]=multi_wiltheta(M,K,C,d,e,f,P,theta,t)
n=length(P(:,1));m=length(P(1,:));
u=zeros(n,m);u1=zeros(n,m);u2=zeros(n,m);uu=zeros(n,m);PP=zeros(n,m);
u(:,1)=d;u1(:,1)=e;u2(:,1)=f;
l=6/(theta*t)^2;q=3/(theta*t);r=6/(theta*t);s=theta*t/2;
KK=K+l*M+q*C;
for i=1:m-1
pp(:,i)=P(:,i)+theta*(P(:,i+1)-P(:,i))+M*(l*u(:,i)+r*u1(:,i)+2*u2(:,i))+C*(q*u(:,i)+2*u1(:,i)+s*u2(:,i));
uu(:,i)=KK\pp(:,i);
u2(:,i+1)=u2(:,i)+l/theta*(uu(:,i)-u(:,i))-l*t*u1(:,i)-3/theta*u2(:,i);
u1(:,i+1)=u1(:,i)+t/2*(u2(:,i+1)+u2(:,i));
u(:,i+1)=u(:,i)+t*u1(:,i)+t^2/6*(u2(:,i+1)+2*u2(:,i));
end
end