FEM一维

 

  

问题原型 

边界条件：

 

（1）取 alpha，beta=1 f=x

(2)p r q=1

(4)(5)

实际结果如下：

 

 完整代码如下

clc; clear all
M=40; % Local-segment number
N=M+1;%points number
av=1;bv=1;L=1;le=L/M;
alpha=ones(M,1)'*av;
beta=ones(M,1)'*bv;
l=ones(M,1)'*le;
f=[1:M]*le;
% 2. boundary condition
p=1;%x=0 position
r=1;%x=L position
q=1;
%3 print input

%4 get K: K(i,i) to a , k(i+1,i) to c
i=2:N-1;
a=[alpha(1)/l(1)+beta(1)*l(1)/3,...
        alpha(i-1)./l(i-1)+beta(i-1).*l(i-1)/3  ...
          + alpha(i)./l(i)+beta(i).*l(i)/3, ...
          alpha(M)/l(M)+beta(M)*l(M)/3]';
i=1:N-1;
c=[-alpha(i)./l(i)+beta(i).*l(i)/6]';

%5 get b
i=2:N-1;
b=[f(1)*l(1)/2,  ...
    f(i-1).*l(i-1)/2+f(i).*l(i)/2, ...
    f(M)*l(M)/2]';

%6 boudary condition apply

a(end)=a(end)+r; %x=L third condition
b(end)=b(end)+q;

K=diag(a,0)+diag(c,-1)+diag(c,+1);
K(1,1)=1; b(1)=p;K(1,2:end)=0;%x=0 1st condition

%7 compute
x=(0:M)*le;
y=K^(-1)*b;
plot(x,y);

hold on;
x1=(0:100)*0.01;
c1=-1/2/exp(1);
c2=1-c1;
y1=c1*exp(x1)+c2*exp(-x1)+x1;
plot(x1,y1,'-r');
hold off;