pardg1cg1.m

various matlab programs to slove various problems

function [u,x]=pardg1cg1(I,n,u0,f,nu,bc)
nx = n(1);          
h = (I(2)-I(1))/nx; 
x  = [I(1):h:I(2)]; 
t = I(3);   
um = (eval(u0))';
nt = n(2);          
k = (I(4)-I(3))/nt;
e = ones(nx+1,1);
Add = spdiags([(h/12-k*nu/(3*h))*e, (h/3+2*k*nu/(3*h))*e, ...
      (h/12-k*nu/(3*h))*e],-1:1,nx+1,nx+1);
Aud = spdiags([(h/12-k*nu/(6*h))*e, (h/3+k*nu/(3*h))*e, ...
      (h/12-k*nu/(6*h))*e],-1:1,nx+1,nx+1);
Ald = spdiags([(-h/12-k*nu/(6*h))*e, (-h/3+k*nu/(3*h))*e, ...
      (-h/12-k*nu/(6*h))*e],-1:1,nx+1,nx+1);
B = spdiags([h*e/6, 2*h*e/3, h*e/6],-1:1,nx+1,nx+1);
Add(1,1) = 1; Add(1,2) = 0; B(1,1)   = 0; B(1,2)   = 0;
Aud(1,1) = 0; Aud(1,2) = 0; Ald(1,1) = 0; Ald(1,2) = 0;
Add(nx+1,nx+1)=1; 
Add(nx+1,nx)=0;
B(nx+1,nx+1)=0; 
B(nx+1,nx) = 0;
Ald(nx+1,nx+1)=0; 
Ald(nx+1,nx)=0;
Aud(nx+1,nx+1)=0; 
Aud(nx+1,nx)=0;
[L,U]=lu([Add Aud; Ald Add]);
x = [I(1)+h:h:I(2)-h]; 
xx=[I(1),x,I(2)];
for time = I(3)+k:k:I(4)
  t = time;     
  fq1 = 0.5*k*h*eval(f);
  t = time-k;   
  fq0 = 0.5*k*h*eval(f);
  rhs0 = [bc(1), fq0, bc(2)];   
  rhs1 = [bc(1), fq1, bc(2)];
  b = [rhs0'; rhs1'] + [B*um; zeros(nx+1,1)];
  y = L \ b;   
  u = U \ y;  
  um = u(nx+2:2*nx+2,1);
end
x = [I(1):h:I(2)]; 
u = um;
return