jacobian.m

fem nonlinear source code using weighted residuals

function [J]=Jacobian(x,c,alpha,N)

J = spalloc(N,N,5*N);
% disp('C in jaco is');
% disp(c)

for i = 1
        %%% small Jacobians %%%
    j(1,1) = 2*(1/(x(i+1)-x(i)))*(0.25)*((5/9)+(8/9)+(5/9))+...
            ((x(i+1)-x(i))/2)*(alpha*0.125)*(2*c(i))*((5/9)*(1+sqrt(3/5))^3+(8/9)+(5/9)*(1-sqrt(3/5))^3)+...
            ((x(i+1)-x(i))/2)*(0.125*alpha)*(2*c(i+1))*((5/9)*(1-sqrt(3/5))*(1+sqrt(3/5))^2+(8/9)+(5/9)*(1+sqrt(3/5))*(1-sqrt(3/5))^2);
            
    j(1,2) = 2*(1/(x(i+1)-x(i)))*(-0.25*((5/9)+(8/9)+(5/9)))+...
            ((x(i+1)-x(i))/2)*(0.125*alpha*2*c(i))*((5/9)*(1-sqrt(3/5))*(1+sqrt(3/5))^2+(8/9)+(5/9)*(1+sqrt(3/5))*(1-sqrt(3/5))^2)+...
            ((x(i+1)-x(i))/2)*(alpha*0.125*2*c(i+1))*((5/9)*(1-sqrt(3/5))^2*(1+sqrt(3/5))+(8/9)+(5/9)*(1+sqrt(3/5))^2*(1-sqrt(3/5)));
    
    j(2,1)= 2*(1/(x(i+1)-x(i)))*(-0.25)*((5/9)+(8/9)+(5/9))+...
            ((x(i+1)-x(i))/2)*(0.125*alpha*2*c(i))*((5/9)*(1-sqrt(3/5))*(1+sqrt(3/5))^2+(8/9)+(5/9)*(1+sqrt(3/5))*(1-sqrt(3/5))^2)+...
            ((x(i+1)-x(i))/2)*(0.125*alpha*2*c(i+1))*((5/9)*(1+sqrt(3/5))*(1-sqrt(3/5))^2+(8/9)+(5/9)*(1-sqrt(3/5))*(1+sqrt(3/5))^2);
        
    j(2,2) = 2*(1/(x(i+1)-x(i)))*(0.25)*((5/9)+(8/9)+(5/9))+...
            ((x(i+1)-x(i))/2)*(0.125*alpha*2*c(i))*((5/9)*(1+sqrt(3/5))*(1-sqrt(3/5))^2+(8/9)+(5/9)*(1-sqrt(3/5))*(1+sqrt(3/5))^2)+...
            ((x(i+1)-x(i))/2)*(alpha*0.125*2*c(i+1))*((5/9)*(1-sqrt(3/5))^3+(8/9)+(5/9)*(1+sqrt(3/5))^3);   
end
disp(j);

    J(1,1) = j(1,1);
    
for k = 1:(N-1)
    J(k,k) = J(k,k)+j(1,1);
    J(k,k+1) = j(1,2);
    J(k+1,k) = j(2,1);
    J(k+1,k+1) = J(k+1,k+1)+j(2,2);
end
    J(1,:) = 0;
    J(1,1)=1;
    
end