% file = irks.m % % This function computes explicit and implicit % irks methods. % % lambda=0 explicit, with error constant cp1 % lambda>0 Lstable, cp1 is ignored % % function glm=irks(lambda,cp1,c,beta,T,dbeta,dT) % % Date: 17 March 2006 % % ----------------------------------------------- function glm=irks(lambda,cp1,c,beta,T,dbeta,dT) s=length(c); p=s-1; if lambda>1e-5 % compute the error constant for L-stabiliy pp=factorial(p)./(factorial(0:p).^2.*factorial(p:-1:0)); cp1=polyval([pp,1/factorial(s)],-lambda); end colx=[0;0;-beta']; K=diag(ones(s-1,1),1); %E=polyvalm(1./gamma(s:-1:1),K); E=toeplitz([1;zeros(s-1,1)],1./gamma(1:s)); Psi=zeros(s,s); Psi(s,s)=1; for i=s-1:-1:1 Psi(i,s)=lambda*Psi(i+1,s)-colx(i+1); end for i=s-1:-1:1 jj=s-1:-1:1; Psi(i,jj)=lambda*Psi(i+1,jj)+Psi(i+1,jj+1); end Psii=inv(Psi); X=lambda*eye(s)+Psi*K'*Psii; C=repmat(c',1,s).^repmat(0:(s-1),s,1)/diag(gamma(1:s)); es=[zeros(s-1,1);1]; Xes=X(:,s); Ees=E(:,s); % Generate the linear systems for Btilde rhs=-[ Psii(2:end,:)*(E*Xes-X*Ees) zeroV(E,T) 1/factorial(s)-cp1 zeros(s-1,1)]; Bt=zeros(s,s); G=[]; for i=1:s for j=1:i Bt(i,j)=1; eq=[ Bt(2:end,:)*C*(es-K*Xes)+Psii(2:end,:)*(E*Xes-X*Ees) % IRKS zeroV(E-Psi*Bt*C*K,T) % ZeroStab 1/factorial(s)-cp1-Bt(1,1)*(C(1,:)-lambda*C(1,:)*K)*Psi(:,s) % Error Const (Psi(2,:)*Bt(:,1:s-1))' % FSAL condition ]; G=[G,eq+rhs]; Bt(i,j)=0; end end bt=G\rhs; for i=1:s l=1+i*(i-1)/2; jj=l:l+i-1; Bt(i,1:i)=bt(jj)'; end glm.btilde=Bt; glm.B=Psi*Bt; glm.A=(Bt\(Psi\X))*glm.B; glm.V=E-Psi*Bt*C*K; glm.U=C-glm.A*C*K; if nargin<=5, return, end % ==================================== dcolx=[0;0;-dbeta']; dPsi=zeros(s,s); for i=s-1:-1:1 dPsi(i,s)=lambda*dPsi(i+1,s)-dcolx(i+1); end for i=s-1:-1:1 jj=s-1:-1:1; dPsi(i,jj)=lambda*dPsi(i+1,jj)+dPsi(i+1,jj+1); end dPsii=-inv(Psi)*dPsi*inv(Psi); dX=dPsi*K'*inv(Psi)+Psi*K'*dPsii; dXes=dX(:,s); row2end=@(r)reshape(r(2:end),s-1,1); dEQ1=@(Bt)row2end(... Bt*C*(-K*dXes)+dPsii*(E*Xes-X*Ees)+Psi\(E*dXes-dX*Ees)); % IRKS condition dEQ2=@(Bt)zerodV(E-Psi*Bt*C*K,-dPsi*Bt*C*K,T,dT); % Zero stability dEQ3=@(Bt)Bt(1,1)*(C(1,:)-lambda*C(1,:)*K)*Psi(:,s); % Error Const dEQ4=@(Bt)(dPsi(2,:)*Bt(:,1:s-1))'; dEQ=@(Bt)[dEQ1(Bt);dEQ2(Bt);dEQ3(Bt);dEQ4(Bt)]; dBt=zeros(s,s); drhs=-dEQ(dBt); dG=[]; for i=1:s for j=1:i dBt(i,j)=1; dG=[dG,dEQ(dBt)+drhs]; dBt(i,j)=0; end end dGi=-inv(G)*dG*inv(G); bt=G\drhs+dGi*rhs; for i=1:s l=1+i*(i-1)/2; jj=l:l+i-1; dBt(i,1:i)=bt(jj); end dbtilde=dBt; dB=dPsi*Bt+Psi*dBt; dBi=-inv(glm.B)*dB*inv(glm.B); dA=dBi*X*glm.B+inv(glm.B)*dX*glm.B+inv(glm.B)*X*dB; dV=-dPsi*Bt*C*K-Psi*dBt*C*K; dU=-dA*C*K; glm.btilde=dbtilde; glm.A=dA; glm.B=dB; glm.U=dU; glm.V=dV; %-------- some helpers ....................... function elem=zeroV(V,T) elem=[trilpart(T\V(3:end,3:end)*T)]; function elem=zerodV(V,dV,T,dT) dTi=-inv(T)*dT*inv(T); elem=[trilpart(dTi*V(3:end,3:end)*T + ... T\dV(3:end,3:end)*T + ... T\V(3:end,3:end)*dT)]; function elem=trilpart(T) s=length(T); elem=zeros(s*(s+1)/2,1); for i=1:s jj=1:i; kk=i*(i-1)/2+jj; elem(kk)=T(i,jj); end