residual.m

来自「fem nonlinear source code using weighted」· M 代码 · 共 35 行

function [R]=Residual(x,c,alpha,N,R,R1,R2)

for i=1
    R1(i) = 2*(1/(x(i+1)-x(i)))*(0.25*(c(i)-c(i+1))*((5/9)+(8/9)+(5/9)))+...
            ((x(i+1)-x(i))/2)*(alpha*0.125*c(i)^2)*((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)*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)+...
            ((x(i+1)-x(i))/2)*(alpha*0.125*c(i+1)^2)*((5/9)*(1-sqrt(3/5))^2*(1+sqrt(3/5))+(8/9)+(5/9)*(1+sqrt(3/5))^2*(1-sqrt(3/5)));
        
        %%% R2 %%%
    R2(i) = 2*(1/(x(i+1)-x(i)))*(0.25*(c(i+1)-c(i))*((5/9)+(8/9)+(5/9)))+...
            ((x(i+1)-x(i))/2)*(0.125*alpha*c(i)^2)*((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)*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)+...
            ((x(i+1)-x(i))/2)*(alpha*0.125*c(i+1)^2)*((5/9)*(1-sqrt(3/5))^3+(8/9)+(5/9)*(1+sqrt(3/5))^3);
    R(i)=R(i)+R1(i); R(i+1)=R(i+1)+R2(i);
end

for i = 2:(N-1) %looping over elements
    
    R1(i) = 2*(1/(x(i+1)-x(i)))*(0.25*(c(i)-c(i+1))*((5/9)+(8/9)+(5/9)))+...
            ((x(i+1)-x(i))/2)*(alpha*0.125*c(i)^2)*((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)*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)+...
            ((x(i+1)-x(i))/2)*(alpha*0.125*c(i+1)^2)*((5/9)*(1-sqrt(3/5))^2*(1+sqrt(3/5))+(8/9)+(5/9)*(1+sqrt(3/5))^2*(1-sqrt(3/5)));
        
        %%% R2 %%%
    R2(i) = 2*(1/(x(i+1)-x(i)))*(0.25*(c(i+1)-c(i))*((5/9)+(8/9)+(5/9)))+...
            ((x(i+1)-x(i))/2)*(0.125*alpha*c(i)^2)*((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)*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)+...
            ((x(i+1)-x(i))/2)*(alpha*0.125*c(i+1)^2)*((5/9)*(1-sqrt(3/5))^3+(8/9)+(5/9)*(1+sqrt(3/5))^3);
        
        %%% Composite residual %%%
    R(i)=R(i)+R1(i); R(i+1)=R(i+1)+R2(i);
  
end
    R(1)=0;
end