youxianyuan2.m

来自「数值分析方法,主要里面涵盖了计算方法方面matlab程序！」· M 代码 · 共 266 行

clear
clc
x=0:1:2;
y=0:0.5:2;
k=0;
N1=length(x);
N2=length(y);
eps=1e-5;
S=1/4;
KK=zeros(15);
N=0;
fprintf('选取节点的方案:\n');
for j=1:length(y)
    for i=1:length(x)
        k=k+1;
        fprintf('第%d个节点坐标:\n',k);
        vj(k,:)=[x(i),y(j)];
        fprintf('(%0.2f,%0.2f)\n',x(i),y(j));
    end
end
for t=1:k
    if vj(t,1)~=0&vj(t,2)~=0&vj(t,1)~=2&vj(t,2)~=2
        N=N+1;
        if rem(N,2)==0
            continue;
        end
        fprintf('第%d个中心节点为%d,相应单元号的节点号为:\n',N,t);
        for i=t-N1-1:t-N1
            fprintf('(%d,%d,%d)\n',t,i,i+1);
            y=[vj(t,2),vj(i,2),vj(i+1,2)];
            a(1)=y(2)-y(3);   
            a(2)=y(3)-y(1);
            a(3)=y(1)-y(2);
            x=[vj(t,1),vj(i,1),vj(i+1,1)];
            b(1)=x(3)-x(2);   
            b(2)=x(1)-x(3);
            b(3)=x(2)-x(1);
            if N==1
                K=[a',b'];
                Ke=1/(4*S)*K*K';
                KK(t,t)=KK(t,t)+Ke(1,1);
                KK(i,i)=KK(i,i)+Ke(2,2);
                KK(i+1,i+1)=KK(i+1,i+1)+Ke(3,3);
                KK(t,i)=KK(t,i)+Ke(1,2);
                KK(i,t)=KK(i,t)+Ke(2,1);
                KK(i+1,i)=KK(i+1,i)+Ke(3,2);
                KK(t,i+1)=KK(t,i+1)+Ke(1,3);
                KK(i+1,t)=KK(i+1,t)+Ke(3,1);
                KK(i,i+1)=KK(i,i+1)+Ke(2,3);
            end
            if N==3  
                K=[a',b'];
                Ke=1/(4*S)*K*K';
                KK(t,t)=KK(t,t)+Ke(1,1);
                KK(i,i)=KK(i,i)+Ke(2,2);
                KK(i+1,i+1)=KK(i+1,i+1)+Ke(3,3);
                KK(t,i)=KK(t,i)+Ke(1,2);
                KK(i,t)=KK(i,t)+Ke(2,1);
                KK(i+1,i)=KK(i+1,i)+Ke(3,2);
                KK(t,i+1)=KK(t,i+1)+Ke(1,3);
                KK(i+1,t)=KK(i+1,t)+Ke(3,1);
                KK(i,i+1)=KK(i,i+1)+Ke(2,3);
            end
        end
        fprintf('(%d,%d,%d)\n',t,t-1,t-N1-1); 
        y=[vj(t,2),vj(t-1,2),vj(t-N1-1,2)];
        a(1)=y(2)-y(3);    
        a(2)=y(3)-y(1);
        a(3)=y(1)-y(2);
        x=[vj(t,1),vj(t-1,1),vj(t-N1-1,1)];
        b(1)=x(3)-x(2);   
        b(2)=x(1)-x(3);
        b(3)=x(2)-x(1);
        if N==1
            K=[a',b'];
            Ke=1/(4*S)*K*K';
            KK(t,t)=KK(t,t)+Ke(1,1);
            KK(t-1,t-1)=KK(t-1,t-1)+Ke(2,2);
            KK(t-N1-1,t-N1-1)=KK(t-N1-1,t-N1-1)+Ke(3,3);
            KK(t,t-1)=KK(t,t-1)+Ke(1,2);
            KK(t-1,t)=KK(t-1,t)+Ke(2,1);
            KK(t-N1-1,t-1)=KK(t-N1-1,t-1)+Ke(3,2);
            KK(t,t-N1-1)=KK(t,t-N1-1)+Ke(1,3);
            KK(t-N1-1,t)=KK(t-N1-1,t)+Ke(3,1);
            KK(t-1,t-N1-1)=KK(t-1,t-N1-1)+Ke(2,3);
        end
        if N==3
            K=[a',b'];
            Ke=1/(4*S)*K*K';
            KK(t,t)=KK(t,t)+Ke(1,1);
            KK(t-1,t-1)=KK(t-1,t-1)+Ke(2,2);
            KK(t-N1-1,t-N1-1)=KK(t-N1-1,t-N1-1)+Ke(3,3);
            KK(t,t-1)=KK(t,t-1)+Ke(1,2);
            KK(t-1,t)=KK(t-1,t)+Ke(2,1);
            KK(t-N1-1,t-1)=KK(t-N1-1,t-1)+Ke(3,2);
            KK(t,t-N1-1)=KK(t,t-N1-1)+Ke(1,3);
            KK(t-N1-1,t)=KK(t-N1-1,t)+Ke(3,1);
            KK(t-1,t-N1-1)=KK(t-1,t-N1-1)+Ke(2,3);
        end       
        fprintf('(%d,%d,%d)\n',t,t+N1-1,t-1); 
        y=[vj(t,2),vj(t+N1-1,2),vj(t-1,2)];
        a(1)=y(2)-y(3);    
        a(2)=y(3)-y(1);
        a(3)=y(1)-y(2);
        x=[vj(t,1),vj(t+N1-1,1),vj(t-1,1)];
        b(1)=x(3)-x(2);   
        b(2)=x(1)-x(3);
        b(3)=x(2)-x(1);
        if N==1
            K=[a',b'];
            Ke=1/(4*S)*K*K';
            KK(t,t)=KK(t,t)+Ke(1,1);
            KK(t+N1-1,t+N1-1)=KK(t+N1-1,t+N1-1)+Ke(2,2);
            KK(t-1,t-1)=KK(t-1,t-1)+Ke(3,3);
            KK(t,t+N1-1)=KK(t,t+N1-1)+Ke(1,2);
            KK(t+N1-1,t)=KK(t+N1-1,t)+Ke(2,1);
            KK(t-1,t+N1-1)=KK(t-1,t+N1-1)+Ke(3,2);
            KK(t,t-1)=KK(t,t-1)+Ke(1,3);
            KK(t-1,t)=KK(t-1,t)+Ke(3,1);
            KK(t+N1-1,t-1)=KK(t+N1-1,t-1)+Ke(2,3);
        end
        if N==3
            K=[a',b'];
            Ke=1/(4*S)*K*K';
            KK(t,t)=KK(t,t)+Ke(1,1);
            KK(t+N1-1,t+N1-1)=KK(t+N1-1,t+N1-1)+Ke(2,2);
            KK(t-1,t-1)=KK(t-1,t-1)+Ke(3,3);
            KK(t,t+N1-1)=KK(t,t+N1-1)+Ke(1,2);
            KK(t+N1-1,t)=KK(t+N1-1,t)+Ke(2,1);
            KK(t-1,t+N1-1)=KK(t-1,t+N1-1)+Ke(3,2);
            KK(t,t-1)=KK(t,t-1)+Ke(1,3);
            KK(t-1,t)=KK(t-1,t)+Ke(3,1);
            KK(t+N1-1,t-1)=KK(t+N1-1,t-1)+Ke(2,3);
        end
        for i=t+N1+1:-1:t+N1
            fprintf('(%d,%d,%d)\n',t,i,i-1);
            y=[vj(t,2),vj(i,2),vj(i-1,2)];
            a(1)=y(2)-y(3);   
            a(2)=y(3)-y(1);
            a(3)=y(1)-y(2);
            x=[vj(t,1),vj(i,1),vj(i-1,1)];
            b(1)=x(3)-x(2);   
            b(2)=x(1)-x(3);
            b(3)=x(2)-x(1);        
            if N==1
                K=[a',b'];
                Ke=1/(4*S)*K*K';
                KK(t,t)=KK(t,t)+Ke(1,1);
                KK(i,i)=KK(i,i)+Ke(2,2);
                KK(i-1,i-1)=KK(i-1,i-1)+Ke(3,3);
                KK(t,i)=KK(t,i)+Ke(1,2);
                KK(i,t)=KK(i,t)+Ke(2,1);
                KK(i-1,i)=KK(i-1,i)+Ke(3,2);
                KK(t,i-1)=KK(t,i-1)+Ke(1,3);
                KK(i-1,t)=KK(i-1,t)+Ke(3,1);
                KK(i,i-1)=KK(i,i-1)+Ke(2,3);               
            end
            if N==3
                K=[a',b'];
                Ke=1/(4*S)*K*K';
                KK(t,t)=KK(t,t)+Ke(1,1);
                KK(i,i)=KK(i,i)+Ke(2,2);
                KK(i-1,i-1)=KK(i-1,i-1)+Ke(3,3);
                KK(t,i)=KK(t,i)+Ke(1,2);
                KK(i,t)=KK(i,t)+Ke(2,1);
                KK(i-1,i)=KK(i-1,i)+Ke(3,2);
                KK(t,i-1)=KK(t,i-1)+Ke(1,3);
                KK(i-1,t)=KK(i-1,t)+Ke(3,1);
                KK(i,i-1)=KK(i,i-1)+Ke(2,3);       
            end
        end
        fprintf('(%d,%d,%d)\n',t,t+1,t+N1+1);  
        y=[vj(t,2),vj(t+1,2),vj(t+N1+1,2)];
        a(1)=y(2)-y(3);   
        a(2)=y(3)-y(1);
        a(3)=y(1)-y(2);
        x=[vj(t,1),vj(t+1,1),vj(t+N1+1,1)];
        b(1)=x(3)-x(2);   
        b(2)=x(1)-x(3);
        b(3)=x(2)-x(1);
        if N==1
            K=[a',b'];
            Ke=1/(4*S)*K*K';
            KK(t,t)=KK(t,t)+Ke(1,1);
            KK(t+1,t+1)=KK(t+1,t+1)+Ke(2,2);
            KK(t+N1+1,t+N1+1)=KK(t+N1+1,t+N1+1)+Ke(3,3);
            KK(t,t+1)=KK(t,t+1)+Ke(1,2);
            KK(t+1,t)=KK(t+1,t)+Ke(2,1);
            KK(t+N1+1,t+1)=KK(t+N1+1,t+1)+Ke(3,2);
            KK(t,t+N1+1)=KK(t,t+N1+1)+Ke(1,3);
            KK(t+N1+1,t)=KK(t+N1+1,t)+Ke(3,1);
            KK(t+1,t+N1+1)=KK(t+1,t+N1+1)+Ke(2,3);   
        end
        if N==3
            K=[a',b'];
            Ke=1/(4*S)*K*K';
            KK(t,t)=KK(t,t)+Ke(1,1);
            KK(t+1,t+1)=KK(t+1,t+1)+Ke(2,2);
            KK(t+N1+1,t+N1+1)=KK(t+N1+1,t+N1+1)+Ke(3,3);
            KK(t,t+1)=KK(t,t+1)+Ke(1,2);
            KK(t+1,t)=KK(t+1,t)+Ke(2,1);
            KK(t+N1+1,t+1)=KK(t+N1+1,t+1)+Ke(3,2);
            KK(t,t+N1+1)=KK(t,t+N1+1)+Ke(1,3);
            KK(t+N1+1,t)=KK(t+N1+1,t)+Ke(3,1);
            KK(t+1,t+N1+1)=KK(t+1,t+N1+1)+Ke(2,3);   
        end
        fprintf('(%d,%d,%d)\n',t,t-N1+1,t+1);  
        y=[vj(t,2),vj(t-N1+1,2),vj(t+1,2)];
        a(1)=y(2)-y(3);   
        a(2)=y(3)-y(1);
        a(3)=y(1)-y(2);
        x=[vj(t,1),vj(t-N1+1,1),vj(t+1,1)];
        b(1)=x(3)-x(2);   
        b(2)=x(1)-x(3);
        b(3)=x(2)-x(1);
        if N==1
            K=[a',b'];
            Ke=1/(4*S)*K*K';
            KK(t,t)=KK(t,t)+Ke(1,1);
            KK(t-N1+1,t-N1+1)=KK(t-N1+1,t-N1+1)+Ke(2,2);
            KK(t+1,t+1)=KK(t+1,t+1)+Ke(3,3);
            KK(t,t-N1+1)=KK(t,t-N1+1)+Ke(1,2);
            KK(t-N1+1,t)=KK(t-N1+1,t)+Ke(2,1);
            KK(t+1,t-N1+1)=KK(t+1,t-N1+1)+Ke(3,2);
            KK(t,t+1)=KK(t,t+1)+Ke(1,3);
            KK(t+1,t)=KK(t+1,t)+Ke(3,1);
            KK(t-N1+1,t+1)=KK(t-N1+1,t+1)+Ke(2,3);
        end
        if N==3
            K=[a',b'];
            Ke=1/(4*S)*K*K';
            KK(t,t)=KK(t,t)+Ke(1,1);
            KK(t-N1+1,t-N1+1)=KK(t-N1+1,t-N1+1)+Ke(2,2);
            KK(t+1,t+1)=KK(t+1,t+1)+Ke(3,3);
            KK(t,t-N1+1)=KK(t,t-N1+1)+Ke(1,2);
            KK(t-N1+1,t)=KK(t-N1+1,t)+Ke(2,1);
            KK(t+1,t-N1+1)=KK(t+1,t-N1+1)+Ke(3,2);
            KK(t,t+1)=KK(t,t+1)+Ke(1,3);
            KK(t+1,t)=KK(t+1,t)+Ke(3,1);
            KK(t-N1+1,t+1)=KK(t-N1+1,t+1)+Ke(2,3);
        end
    end
end    
u(1)=50;
u(2)=50;
u(3)=50;
u(13)=100;
u(14)=100;
u(15)=100;
u;
K11=KK(1:N1,1:N1);
K12=KK(1:N1,N1+1:15-N1);
K13=KK(1:N1,16-N1:15);   
K21=KK(N1+1:15-N1,1:N1);
K22=KK(N1+1:15-N1,N1+1:15-N1);
K23=KK(N1+1:15-N1,16-N1:15);
K31=KK(16-N1:15,1:N1);   
K32=KK(16-N1:15,N1+1:15-N1);
K33=KK(16-N1:15,16-N1:15);
Y=-K21*u(1:3)'-K23*u(13:15)';
X=inv(K22)*Y;
fprintf('边界结果为：\n')
n=size(K22,1);
for i=4:n+3
    fprintf('u[%d]=%f\n',i,X(i-3));
end