fem.m

1、 写出此两点边值问题的积分形式； 2、 构造有限元方程； 3、 编程计算显示计算结果。

%求解一维问题的有限元法的通用子程序
%输入变量解释：a,b为变量的初末值，n为等分的份数，mo,mn为函数的初末值
function u=fem(a,b,n,m0,mn)
h=(b-a)/n;
%用MATLAB的内部函数来求解有限元方程系数的积分，b的求解是用外部定义的子函数求解
syms h1 kesi1;
f1=1/h+h*(1-kesi1)^2;
f1=int(f1,kesi1,0,1);
ai1i1=eval(f1);
syms h2 kesi2;
f2=-1/h+h*kesi2*(1-kesi2);
f2=int(f2,kesi2,0,1);
h2=h;
ai1i=eval(f2);
syms h3 kesi3;
f3=1/h+h*kesi3^2;
f3=int(f3,kesi3,0,1);
h3=h;
aii=eval(f3);
%求解有限元方程的系数矩阵和右端的常数向量
A=zeros(n+1,n+1);
B=zeros(n+1);
for i=2:n+1
    A1=zeros(n+1,n+1);
    B1=zeros(n+1);
    A1(i-1,i-1)=ai1i1;
    A1(i-1,i)=ai1i;
    A1(i,i-1)=ai1i;
    A1(i,i)=aii;
    A=A+A1;
    x(i-1)=h*(i-2);
    B1(i-1)=b1(x(i-1),h);
    B1(i)=b2(x(i-1),h);
    B=B+B1;
end
%调整有限元方程的系数矩阵和右端的常数向量
A=A(2:n,2:n);
B=B(2:n);
x=x(2:n);
%用不选主元的gauss消去法求解有限元方程
[L,U]=nongauss(A);
y=downsan(L,B');
u=upsan(U,y);
%把函数值变为变换前的函数值
u=u+x'*(exp(1)-exp(-1));