vgpjac

 jacobian for vector GP

0001 function [c,fu,flam,b]=vgpjac(p,u,lam) 
0002 % jacobian for vector GP
0003 mu1=p.mu1;mu2=p.mu2; c=[1;0;0;1;1;0;0;1;1;0;0;1;1;0;0;1]; %mu=lam;
0004 flam=[0;0;0;0]; % do numerically since inexpensive
0005 x=p.points(1,:)'; y=p.points(2,:)'; r=x.^2+y.^2; 
0006 pot=p.pa*pdeintrp(p.points,p.tria,r);
0007 ui=pdeintrp(p.points,p.tria,u);
0008 u1=ui(1,:); v1=ui(2,:);u2=ui(3,:); v2=ui(4,:); 
0009 a1=u1.^2+v1.^2; a2=u2.^2+v2.^2; 
0010 a1u1=2*u1; a1v1=2*v1; a2u2=2*u2; a2v2=2*v2; 
0011 f1u1=mu1-pot+a1u1.*u1+a1+0.5*a2; 
0012 f1v1=a1v1.*u1; f1u2=0.5*a2u2.*u1; f1v2=0.5*a2v2.*u1; 
0013 f2v1=mu1-pot+a1v1.*v1+a1+0.5*a2; 
0014 f2u1=a1u1.*v1; f2u2=0.5*a2u2.*v1; f2v2=0.5*a2v2.*v1; 
0015 f3u2=mu2-pot+a2u2.*u2+a2+0.5*a1;
0016 f3u1=0.5*a1u1.*u2; f3v1=0.5*a1v1.*u2; f3v2=a2v2.*u2; 
0017 f4v2=mu2-pot+a2v2.*v2+a2+0.5*a1;
0018 f4u1=0.5*a1u1.*v2; f4v1=0.5*a1v1.*v2; f4u2=a2u2.*v2;
0019 fu=[f1u1; f2u1; f3u1; f4u1; f1v1; f2v1; f3v1; f4v1; ...
0020     f1u2; f2u2; f3u2; f4u2; f1v2; f2v2; f3v2; f4v2]; 
0021 xi=pdeintrp(p.points,p.tria,x);yi=pdeintrp(p.points,p.tria,y);
0022 b=zeros(p.neq*p.neq*2,p.nt); 
0023 b(3,:)=yi;b(4,:)=-xi;
0024 b(9,:)=-yi;b(10,:)=xi;
0025 b(23,:)=yi;b(24,:)=-xi;
0026 b(29,:)=-yi;b(30,:)=xi;
0027 b=lam*b; 
0028