00001 function A=myJac(x,n)
00002
00003 global lambda
00004 global re
00005 global unk
00006
00007 A=sparse(n*n*unk,n*n*unk); % corners are 1
00008
00009 h=1/(n-1.0);
00010
00011
00012
00013 for i=1:n*n*unk
00014 A(i,i)=1;
00015 end
00016
00017
00018
00019
00020 j=1;
00021 for i=2:n-1
00022 index=((j-1)*n+i-1)*unk+1;
00023 A(index,index)=1.0;
00024 A(index+1,index+1)=1.0;
00025 A(index+1,index)=-2.0/h^2;
00026 A(index+1,index+unk*n)=2.0/h^2;
00027 end
00028
00029 j=n;
00030 for i=2:n-1
00031 index=((j-1)*n+i-1)*unk+1;
00032 A(index,index)=1.0;
00033 A(index+1,index+1)=1.0;
00034 A(index+1,index)=-2.0/h^2;
00035 A(index+1,index-unk*n)=2.0/h^2;
00036 end
00037
00038 i=1;
00039 for j=2:n-1
00040 index=((j-1)*n+i-1)*unk+1;
00041 A(index,index)=1.0;
00042 A(index+1,index+1)=1.0;
00043 A(index+1,index)=-2.0/h^2;
00044 A(index+1,index+unk)=2.0/h^2;
00045 end
00046
00047 i=n;
00048 for j=2:n-1
00049 index=((j-1)*n+i-1)*unk+1;
00050 A(index,index)=1.0;
00051 A(index+1,index+1)=1.0;
00052 A(index+1,index)=-2.0/h^2;
00053 A(index+1,index-unk)=2.0/h^2;
00054 end
00055
00056 for j=2:n-1
00057
00058 for i=2:n-1
00059
00060 index=((j-1)*n+i-1)*unk+1;
00061
00062 sc=x(index);
00063 sw=x(index-unk);
00064 se=x(index+unk);
00065 sn=x(index+n*unk);
00066 ss=x(index-n*unk);
00067 vc=x(index+1);
00068 vw=x(index+1-unk);
00069 ve=x(index+1+unk);
00070 vn=x(index+1+n*unk);
00071 vs=x(index+1-n*unk);
00072
00073 dSTFdX=(se-sw)/(2.0*h);
00074
00075 dSTFdY=(sn-ss)/(2.0*h);
00076
00077 dVORdX=(ve-vw)/(2.0*h);
00078
00079 dVORdY=(vn-vs)/(2.0*h);
00080
00081
00082
00083 A(index,index)=-4/h^2;
00084 A(index,index-unk)=1/h^2;
00085 A(index,index+unk)=1/h^2;
00086 A(index,index-n*unk)=1/h^2;
00087 A(index,index+n*unk)=1/h^2;
00088 A(index,index+1)=1;
00089
00090 A(index+1,index+1)=4/h^2;
00091 A(index+1,index+1+unk)=re*dSTFdY/(2.0*h)-1/h^2;
00092 A(index+1,index+1-unk)=-re*dSTFdY/(2.0*h)-1/h^2;
00093 A(index+1,index+1+unk*n)=-re*dSTFdX/(2.0*h)-1/h^2;
00094 A(index+1,index+1-unk*n)=re*dSTFdX/(2.0*h)-1/h^2;
00095 A(index+1,index+unk)=-re*dVORdY/(2.0*h);
00096 A(index+1,index-unk)=re*dVORdY/(2.0*h);
00097 A(index+1,index+unk*n)=re*dVORdX/(2.0*h);
00098 A(index+1,index-unk*n)=-re*dVORdX/(2.0*h);
00099
00100 end
00101
00102 end
00103
00104
00105 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00106 % for j=2:n-1
00107 %
00108 % for i=2:n-1
00109 %
00110 % index=(j-1)*n+i;
00111 %
00112 % A(index,index)=-4/h^2+lambda*exp(x(index));
00113 % A(index,index-1)=1/h^2;
00114 % A(index,index+1)=1/h^2;
00115 % A(index,index-n)=1/h^2;
00116 % A(index,index+n)=1/h^2;
00117 %
00118 %
00119 % end
00120 %
00121 % end
