Add files via upload

Author: sandeep026

BVP.m

@@ -0,0 +1,57 @@
+function [err]=BVP(N)
+
+% solve boundary value problems using PS methods
+% ode: x''-tx'-x=0
+% BCs: x(-1)=x(1)=1
+%sol : exp((t^2-1)/2)
+
+import casadi.*
+
+%N=10;
+Nx=1;
+t0=-1;
+tf=1;
+
+opti=casadi.Opti();
+
+X=opti.variable(Nx,N+1);
+
+
+%LGL grid points and quadrature weights
+[tau,wi]=legslb(N+1);
+tau=tau';
+wi=wi';
+% differentiation matrix
+D=legslbdiff(N+1,tau);
+
+
+t=(tf-t0)/2*tau+(tf+t0)/2;
+
+Xd=2/(tf-t0)*(D*X')';
+Xdd=2/(tf-t0)*(D*Xd')';
+
+res=Xdd-t.*Xd-X;
+
+opti.minimize(0)
+opti.subject_to(res(2:end-1)==0)
+opti.subject_to(X(1)==1)
+opti.subject_to(X(N+1)==1)
+
+opti.set_initial(X,zeros(Nx,N+1))
+
+opti.solver('ipopt')
+sol=opti.solve();
+
+
+Xs=sol.value(X);
+Xa=exp((t.^2-1)/2);
+%plot(t,Xs,'d',t,Xa);
+%xlabel('time [s]')
+%ylabel('$x(t)$')
+%figure
+%semilogx((abs(Xs-Xa)));
+%xlabel('time [s]')
+%ylabel('Error')
+
+err=norm(Xs-Xa,2,'fro')
+end

IVP.m

@@ -0,0 +1,52 @@
+% solve boundary value problems using PS methods
+% ode: x''-tx'-x=0
+%sol : exp((t^2-1)/2)
+
+import casadi.*
+
+N=100;
+Nx=1;
+t0=-1;
+tf=1;
+
+opti=casadi.Opti();
+
+X=opti.variable(Nx,N+1);
+
+
+%LGL grid points and quadrature weights
+[tau,wi]=legslb(N+1);
+tau=tau';
+wi=wi';
+% differentiation matrix
+D=legslbdiff(N+1,tau);
+
+
+t=(tf-t0)/2*tau+(tf+t0)/2;
+
+Xd=2/(tf-t0)*(D*X')';
+Xdd=2/(tf-t0)*(D*Xd')';
+
+res=Xdd-t.*Xd-X;
+
+opti.minimize(0)
+opti.subject_to(res(2:end-1)==0)
+opti.subject_to(X(1)==1)
+opti.subject_to(X(N+1)==1)
+
+opti.set_initial(X,zeros(Nx,N+1))
+
+opti.solver('ipopt')
+sol=opti.solve();
+
+
+Xs=sol.value(X);
+Xa=exp((t.^2-1)/2);
+plot(t,Xs,'d',t,Xa);
+xlabel('time [s]')
+ylabel('$x(t)$')
+figure
+semilogx((abs(Xs-Xa)));
+xlabel('time [s]')
+ylabel('Error')
+

README.md

@@ -0,0 +1,45 @@
+# Solving boundary value problems using pseudospectral method
+
+## Requirements
+
+## Problem description
+
+Consider the following boundary value problem from [^1] in equation 1 and its boundary conditione are provided in equation 2.
+
+$$ 
+\begin{equation}
+\ddot{x}-t\dot{x}-x=0
+\end{equation}
+$$
+
+$$ 
+\begin{equation}
+x(t=-1)=x(t=1)=1
+\end{equation}
+$$
+
+
+## Analytical solution
+
+The system has the analytical solution $$ x(t)=\exp{\frac{t^2-1)}{2}} $$.
+
+## Numerical solution
+
+In order to solve the system  numerically, the problem is discretized using pseudospectral method. $$x(t)$$ is approximated using a global polynomial whose values are defined at the LGL collocation nodes. This polynomial is differentiated to the required degree and constraints are imposed. These include the boundary conditions and the satisfaction of the dynamics at the collocation points. Satifying the BC and dynamics at -1 and 1 leads to a overdetermined system. Hence, the BC are only imposed at the end points and dynamics constraints are ignored. This leads to a system of equation with equal number of variables and equations. This can lead a linear/nonlinear system of equation which can be solved.
+
+## Implementation
+
+In the code, the boundary conditions and  the dynamics are treated as nonlinear constraints by solver "IPOPT" and the objective it set to 0.
+
+## Results
+
+The numerical and analytical solution are identical as shown in the figure below. In addition, the spectral convergence of the method is also additional verified by increasing the grid size progressively.
+
+
+## References
+
+[^1] Wang, L. L., Samson, M. D., & Zhao, X. (2014). A well-conditioned collocation method using a pseudospectral integration matrix. SIAM Journal on Scientific Computing, 36(3), A907-A929.
+
+
+
+

convergence.m

@@ -0,0 +1,12 @@
+clear
+
+N=5:1:20;
+E=zeros(1,length(N));
+
+for i=1:length(N)
+E(i)=BVP(N(i));
+end
+
+semilogy(N,E)
+xlabel('grid size [N]')
+ylabel('discretization error')

legslb.m

@@ -0,0 +1,32 @@
+ function [varargout]=legslb(n)
+
+%  x=legslb(n) returns n Legendre-Gauss-Lobatto points with x(1)=-1, x(n)=1
+%  [x,w]= legslb(n) returns n Legendre-Gauss-Lobatto points and weights
+%  Newton iteration method is used for computing nodes
+%  See Page 99 of the book: J. Shen, T. Tang and L. Wang, Spectral Methods:
+%  Algorithms, Analysis and Applications, Springer Series in Compuational
+%  Mathematics, 41, Springer, 2011. 
+%  Use the function: lepoly() 
+%  Last modified on August 30, 2011
+
+ % Compute the initial guess of the interior LGL points
+  nn=n-1; thetak=(4*[1:nn]-1)*pi/(4*nn+2);
+  sigmak=-(1-(nn-1)/(8*nn^3)-(39-28./sin(thetak).^2)/(384*nn^4)).*cos(thetak);
+  ze=(sigmak(1:nn-1)+sigmak(2:nn))/2;  
+  ep=eps*10;                            % error tolerance for stopping iteration
+  ze1=ze+ep+1;
+ 
+ while max(abs(ze1-ze))>=ep,            % Newton's iteration procedure
+      ze1=ze;
+      [dy,y]=lepoly(nn,ze);
+      ze=ze-(1-ze.*ze).*dy./(2*ze.*dy-nn*(nn+1)*y);  % see Page 99 of the book
+ end;                                   % around 6 iterations are required for n=100
+    varargout{1}=[-1,ze,1]';
+ if nargout==1, return; end;
+   
+ % Use the weight expression (3.188) to compute the weights
+  varargout{2}=[2/(nn*(nn+1)),2./(nn*(nn+1)*y.^2),2/(nn*(nn+1))]'; 
+
+
+
+

legslbdiff.m

@@ -0,0 +1,20 @@
+% D=legslbdiff(n,x) returns the first-order differentiation matrix of size
+% n by n, associated with the Legendre-Gauss-Lobatto points x, which may be computed by 
+% x=legslb(n) or x=legslbndm(n). Note: x(1)=-1 and x(n)=1.
+% See Page 110 of the book: J. Shen, T. Tang and L. Wang, Spectral Methods:
+%  Algorithms, Analysis and Applications, Springer Series in Compuational
+%  Mathematics, 41, Springer, 2011. 
+%  Use the function: lepoly() 
+% Last modified on August 31, 2011
+
+function D=legslbdiff(n,x)
+if n==0, D=[]; return; end;
+xx=x;y=lepoly(n-1,xx); nx=size(x); 
+if nx(2)>nx(1), y=y'; xx=x'; end;  %% y is a column vector of L_{n-1}(x_k)
+  D=(xx./y)*y'-(1./y)*(xx.*y)';  %% compute L_{n-1}(x_j) (x_k-x_j)/L_{n-1}(x_k);     
+                                 % 1/d_{kj} for k not= j (see (3.203)) 
+  D=D+eye(n);                    % add the identity matrix so that 1./D can be operated                                     
+  D=1./D; 
+  D=D-eye(n); D(1,1)=-n*(n-1)/4; D(n,n)=-D(1,1);  %update the diagonal entries  
+  return; 
+ 

lepoly.m

@@ -0,0 +1,39 @@
+function [varargout]=lepoly(n,x)
+  
+% lepoly  Legendre polynomial of degree n
+    % y=lepoly(n,x) is the Legendre polynomial
+    % The degree should be a nonnegative integer 
+    % The argument x should be on the closed interval [-1,1]; 
+    % [dy,y]=lepoly(n,x) also returns the values of 1st-order 
+    %  derivative of the Legendre polynomial stored in dy
+% Last modified on August 30, 2011    
+% Verified with the chart in http://keisan.casio.com/has10/SpecExec.cgi
+
+if nargout==1,
+     if n==0, varargout{1}=ones(size(x));  return; end;
+     if n==1, varargout{1}=x; return; end;
+     polylst=ones(size(x)); poly=x;   % L_0(x)=1, L_1(x)=x
+     for k=2:n,                      % Three-term recurrence relation:  
+	   polyn=((2*k-1)*x.*poly-(k-1)*polylst)/k; % kL_k(x)=(2k-1)xL_{k-1}(x)-(k-1)L_{k-2}(x)
+       polylst=poly; poly=polyn;	
+     end;
+     varargout{1}=polyn;
+end;
+
+if nargout==2,
+     if n==0, varargout{2}=ones(size(x)); varargout{1}=zeros(size(x)); return;end;
+     if n==1, varargout{2}=x; varargout{1}=ones(size(x)); return; end;
+
+     polylst=ones(size(x)); pderlst=zeros(size(x));poly=x; pder=ones(size(x));
+     % L_0=1, L_0'=0, L_1=x, L_1'=1
+    for k=2:n,                          % Three-term recurrence relation:  
+      polyn=((2*k-1)*x.*poly-(k-1)*polylst)/k; % kL_k(x)=(2k-1)xL_{k-1}(x)-(k-1)L_{k-2}(x)
+      pdern=pderlst+(2*k-1)*poly; % L_k'(x)=L_{k-2}'(x)+(2k-1)L_{k-1}(x)
+ 	  polylst=poly; poly=polyn;
+	  pderlst=pder; pder=pdern;
+    end;
+      varargout{2}=polyn; varargout{1}=pdern;
+end;
+
+return
+      

lepolym.m

@@ -0,0 +1,44 @@
+function [varargout]=lepolym(n,x)
+  
+% lepolym  Legendre polynomials of degree up to n, i.e., 0,1, ...,n
+    % y=lepolym(n,x) returns the Legendre polynomials
+    % The degree should be a nonnegative integer 
+    % The argument x is a vector be on the closed interval [-1,1]; 
+    % [dy,y]=lepolym(n,x) also returns the values of 1st-order 
+    %  derivatives of the Legendre polynomials upto n stored in dy
+    % Note: y (and likewise for dy) saves L_0(x), L_1(x), ...., L_n(x) by rows
+    % i.e., L_k(x) is the (k+1)th row of the matrix y (or dy)
+% Last modified on August 30, 2011    
+      
+dim=size(x); xx=x; if dim(1)>dim(2), xx=xx'; end; % xx is a row-vector
+if nargout==1,
+     if n==0, varargout{1}=ones(size(xx));  return; end;
+     if n==1, varargout{1}=[ones(size(xx));xx]; return; end;
+     polylst=ones(size(xx)); poly=xx;   % L_0(x)=1, L_1(x)=x
+     y=[polylst;poly];
+     for  k=2:n,                      % Three-term recurrence relation:  
+	   polyn=((2*k-1)*xx.*poly-(k-1)*polylst)/k; % kL_k(x)=(2k-1)xL_{k-1}(x)-(k-1)L_{k-2}(x)
+       polylst=poly; poly=polyn; y=[y;polyn];	
+     end;
+     varargout{1}=y;
+end;
+
+if nargout==2,
+     if n==0, varargout{2}=ones(size(xx)); varargout{1}=zeros(size(xx)); return;end;
+     if n==1, varargout{2}=[ones(size(xx));xx]; 
+              varargout{1}=[zeros(size(xx));ones(size(xx))]; return; end;
+
+     polylst=ones(size(xx)); pderlst=zeros(size(xx));poly=xx; pder=ones(size(xx));
+     y=[polylst;poly]; dy=[pderlst;pder];
+     % L_0=1, L_0'=0, L_1=x, L_1'=1
+    for k=2:n,                          % Three-term recurrence relation:  
+      polyn=((2*k-1)*xx.*poly-(k-1)*polylst)/k; % kL_k(x)=(2k-1)xL_{k-1}(x)-(k-1)L_{k-2}(x)
+      pdern=pderlst+(2*k-1)*poly;  % L_k'(x)=L_{k-2}'(x)+(2k-1)L_{k-1}(x)
+ 	  polylst=poly; poly=polyn; y=[y;polyn];
+	  pderlst=pder; pder=pdern; dy=[dy;pdern];
+    end;
+      varargout{2}=y; varargout{1}=dy;
+end;
+
+return
+