function sol=bvp4cfe(odefun,bcfun,solinit,options) ## Integrates an ODE using the finite element method ## ## sol=bvp4cfe(odefun,bcfun,solinit,options) ## ## Input: ## ## odefun: A function handle of the ODE system, each row corresponds to a first derivative of the function you want to integrate, each column corresponds to a different point. The function must accept vectorized input. ## bcfun: A function handle which accepts two vectors of the first and last points and returns a vector which will equal zero when the boundary conditions are satisfied ## solinit: a structure with the following fields: ## * x: the initial mesh a row vector of ordered points from the beginning to the end of the domain ## * y: the guess of the solution a matrix where each row is a different variable and each column is a different point ## options: a structure with the following fields (optional): ## * RelTol: Relative tolerance, default 1e-3 ## * AbsTol: Absolute tolerance, default 1e-6 ## ## Output: ## ## sol: a structure with the following fields ## * x: the locations the functions are evaluated at ## * y: the values of the functions corresponding to x, each row corresponds to a different function ## * yp: the derivative of y with respect to x formatted the same way as y ## ## Example: ## ## f=@(t,y)[y(2,:);-y(1,:)-.5*y(2,:)]; ## bcfun=@(ya,yb)[ya(1)-1;yb(1)]; ## solinit.x=linspace(0,10,5); ## solinit.y=zeros(2,5); ## sol=bvp4cfe(f,bcfun,solinit); ## plot(sol.x,sol.y); if (isfield(solinit,"x")) if columns(solinit.x)<2 error("bvp4c: initial mesh solinit.x must have 2 or more points."); endif t = solinit.x; else error("bvp4c: missing initial mesh solinit.x"); end if (isfield(solinit,"y")) y = solinit.y; else error("bvp4c: missing initial guess"); end if (isfield(solinit,"parameters")) error("bvp4c: solving for unknown parameters is not yet supported"); end RelTol = 1e-3; AbsTol = 1e-6; if ( nargin > 3 ) if (isfield(options,"RelTol")) RelTol = options.RelTol; endif if (isfield(options,"RelTol")) AbsTol = options.AbsTol; endif endif ##profile off; ##profile clear; ##tic; ##profile on; nip=5;#number of points to use when intergrating [w01,t01]=gausslegendrewx(nip); nvar=size(y,1); nnode=length(t); ydy=zeros(nvar,2*nnode); ydy(:,1:2:end)=y; if (isfield(solinit,"yp")) ydy(:,2:2:end) = solinit.yp; else ydy(:,2:2:end) = odefun(t,y); end ydy=ydy(:); maxref=14; for kk=1:maxref h=diff(t); spt=logical(kron(speye(nnode-1),ones(nip,1))); wint=sparse((nnode-1)*nip,nnode-1); wint(spt)=kron(ones(nnode-1,1),w01).*kron(h.',ones(nip,1)); tint=kron(ones(1,nnode-1),t01.').*kron(h,ones(1,nip))+kron(t(1:nnode-1),ones(1,nip)); telem=kron(ones(1,nnode-1),t01.').*kron(h,ones(1,nip)); t3210=telem.^((3:-1:0).'); t2100=((3:-1:0).').*telem.^([(2:-1:0).';0]); spt=logical(kron(speye(nnode-1),ones(4,nip))); t3210sp=sparse((nnode-1)*4,nip*(nnode-1)); t3210sp(spt)=t3210; t2100sp=sparse((nnode-1)*4,nip*(nnode-1)); t2100sp(spt)=t2100; m=elementinterp(h); phi=m*t3210sp; dphi=m*t2100sp; [ydy,r,ex]=fsolve(@(ydy)resasem(ydy,nvar,nnode,nip,phi,dphi,tint,wint,odefun,bcfun),ydy); ##ex y=reshape(ydy,nvar,[]); y(:,2:2:end)=[]; dy=reshape(ydy,nvar,[]); dy(:,1:2:end)=[]; tm=t(1:end-1)+h/2; [ym dy_est]=evaly(.5,1,nnode,nvar,h,ydy); f_est=odefun(tm,ym); err=max(abs(f_est-dy_est)); AbsErr=max(err); RelErr=AbsErr/norm(dy,inf); if ((AbsErr>=AbsTol)&&(RelErr>=RelTol))&&kk=AbsTol)&(err./max(max(abs(dy)))>=RelTol); t_add=tm(ref_int); y_add=ym(:,ref_int); dy_add=dy_est(:,ref_int); [t ti]=sort([t,t_add]); y=[y y_add](:,ti); dy=[dy dy_add](:,ti); nnode=length(t); ydy=zeros(nvar,2*nnode); ydy(:,1:2:end)=y; ydy(:,2:2:end) = dy; ydy=ydy(:); else break end end sol.solver="bvp4c"; sol.x=t; sol.y=y; sol.yp=dy; ##profile off; ##toc ##profshow () ##profexplore () end function [y dy]=evaly(te,np,nnode,nvar,h,ydy) telem=kron(ones(1,nnode-1),te).*kron(h,ones(1,np)); t3210=telem.^((3:-1:0).'); spt=logical(kron(speye(nnode-1),ones(4,np))); t3210sp=sparse((nnode-1)*4,np*(nnode-1)); t3210sp(spt)=t3210; m=elementinterp(h); phi=m*t3210sp; t2100=((3:-1:0).').*telem.^([(2:-1:0).';0]); t2100sp=sparse((nnode-1)*4,np*(nnode-1)); t2100sp(spt)=t2100; dphi=m*t2100sp; nv=reshape(ydy,nvar,2,nnode); nv=[nv(:,:,1:end-1) nv(:,:,2:end)]; nv=reshape(vec(nv,2),nvar,[]); spt=logical(kron(speye(nnode-1),ones(4,np))); yphi=sparse((nnode-1)*4,np*(nnode-1)); yphi(spt)=phi; y=nv*yphi; ydphi=sparse((nnode-1)*4,np*(nnode-1)); ydphi(spt)=dphi; dy=nv*ydphi; end function m=elementinterp(h) h=permute(h,[1 3 2]); sizeh=size(h); if length(sizeh)==2 sizeh=[sizeh 1]; end zero=zeros(sizeh); one=ones(sizeh); m=[2./h.^3 -3./h.^2 zero one; 1./h.^2 -2./h one zero; -2./h.^3 3./h.^2 zero zero; 1./h.^2 -1./h zero zero]; m=reshape(m,4,4*sizeh(3)); end function res=resasem(ydy,nvar,nnode,nip,phi,dphi,tint,wint,dydt,bcfun) res=[zeros(nnode*2*nvar,1);bcfun(ydy(1:nvar),ydy((nnode-1)*2*nvar+(1:nvar)))]; nv=reshape(ydy,nvar,2,nnode); nv=[nv(:,:,1:end-1) nv(:,:,2:end)]; nv=reshape(vec(nv,2),nvar,[]); spt=logical(kron(speye(nnode-1),ones(4,nip))); yphi=sparse((nnode-1)*4,nip*(nnode-1)); yphi(spt)=phi; ydphi=sparse((nnode-1)*4,nip*(nnode-1)); ydphi(spt)=dphi; ye=nv*yphi; dye=nv*ydphi; r=((kron(ones(4,1),dydt(tint,ye)-dye).*kron(phi,ones(nvar,1)))*wint); res(1:nvar*2)=r(1:nvar*2); adi=(nvar*2+(1:nvar*2*(nnode-2)*2)); le=logical(kron(ones(nnode-2,1),kron([1;0],ones(nvar*2,1)))); re=logical(kron(ones(nnode-2,1),kron([0;1],ones(nvar*2,1)))); res(nvar*2+(1:2*nvar*(nnode-2)))=r(adi(le))+r(adi(re)); res(2*nvar*(nnode-1)+(1:2*nvar))=r(nvar*2+nvar*2*(nnode-2)*2+(1:nvar*2)); end function [w,x]=gausslegendrewx(n) #Calculates the points and weights of Gauss Legendre quadratures with n points for integration between 0 and 1. x=sort(roots(flip(2^n*(bincoeff(n,0:n).*bincoeff((n-1+(0:n))/2,n))))); mp1to1=(n+1:-1:1).'; intval=((1.^mp1to1-(-1).^mp1to1)./mp1to1); w=(x.^(n:-1:0)).'\intval; x=(x+1)/2; w/=2; end %!demo %! t=linspace(0,4,5); %! ic=bvpinit(t,[0;0]); %! dydt=@(t,y)[y(2,:);-abs(y(1,:))]; %! bcfun=@(ya,yb)[ya(1);yb(1)+2]; %! s=bvp4cfe(dydt,bcfun,ic); %! [tt,yy]=ode45(dydt,[0 4],s.y(:,1)); %! plot(s.x,s.y,'-d',tt,yy,'-x'); %!demo %! t=linspace(1/3/pi,1,5); %! ic=bvpinit(t,[0;0]); %! dydt=@(t,y)[y(2,:);-y(1,:)./t.^4-2./t.*y(2,:)]; %! bcfun=@(ya,yb)[ya(1);yb(1)-sin(1)]; %! s=bvp4cfe(dydt,bcfun,ic); %! tt=s.x; %! yy=[sin(1./tt);-cos(1./tt)./tt.^2]; %! disp("maximum error"); %! disp(max(abs(s.y-yy),[],2)); %! [ttt,yyy]=ode45(dydt,[1 1/(3*pi)],s.y(:,end)); %! plot(s.x,s.y,'-d',tt,yy,'-',ttt,yyy,'-o');