线性分组码编程实现_有限元学习笔记9-一维线性元的编程实现(1)

现在来整理一下最终要求解的矩阵还有整个的编程思路。同样直接放上图片。

情形一：齐次初边值条件：

遇到的困难一：本来想用匿名函数直接把函数作为参数传进去来调用p和q，但是我先测试了一下发现不允许这样子。所以现在是要针对具体的问题，先把p和q以及f的函数表达式推导好再直接用。（准备明天再优化一下这个地方，因为针对不同的pq和f都是直接改脚本函数会比较方便）

以该问题为例：

精确解为

推导的原问题的形式为

所以现在是相当于

请读者自行算一下pqf在代入不同的自变量的时候的取值。计算总刚度矩阵时需要用到这个积分。

遇到的困难二：在给定的上面的例子中，程序运行的时候出现了很多问题，第一个问题是得到数值解和给定的精确解函数得到的结果相差甚远，可以说是完全就是求解错了，检查了一遍总刚度矩阵三对角线的构成和代码之后发现没有问题，于是又去检查总载荷向量，最终发现是自己的载荷向量全部都少乘了

和

。

遇到的困难三：把上面的错误改好之后又出现了问题，就是我的第二点的误差非常大。如下图所示：左边是步长为10，右边是步长为100。

在师兄的指导下，最终尝试了用Ab的形式，最后得到的结果看起来像是比较正确的数值解了，但数值解肯定是存在误差的，只能以误差的大小来衡量其数值解的好坏。如下图，同样左边是步长为10的，右边是步长为100的。所以有可能是我才用的专门用来解三对角线方程组的程序有错误。

代码如下：

%得到总刚度矩阵Total stiffness matrix的三条对角线上的元素；
%a为-1对角线上的元素；
%b为主对角线上的元素；
%c为1对角线上的元素；
%h为步长；
function [a,b,c]=TSM(x,a,b,c,h,N);
format long
%计算对角线上的元素值；
for i=1:N-1
    %先计算主对角线上的元素；
    %被积函数表示：p(x(i)+h*s)/h+q(x(i)+h*s)*h*s^2 + p(x(i+1)+h*s)/h+q(x(i+1)+h*s)*h*(1-s)^2 ;
    %还需要额外再算最后一项的积分；
    %func=@s (funp(x(i)+h*s)/h+funq(x(i)+h*s)*h*s^2 + funp(x(i+1)+h*s)/h+funq(x(i+1)+h*s)*h*(1-s)^2);
    %b(i)=int(func,0,1);
    func=@(s) 2/h+((pi^2)*h*(s.^2+(1-s).^2))/4;
    b(i)=quad8(func,0,1);
    %再计算-1对角线上的元素：
    %被积分函数表示：q(x(i+1)+h*s)*h*s*(1-s)-p(x(i+1)+h*s)/h;
    %func=@s (funq(x(i+1)+h*s)*h*s*(1-s)-funp(x(i+1)+h*s)/h);
    %a(i)=int(func,0,1);
    func=@(s) (pi^2*h*s.*(1-s))/4-1/h;
    a(i)=quad8(func,0,1);
end
%func=@s (funp(x(N)+h*s)/h+funq(x(N)+h*s)*h*s^2;
%b(N)=int(func,0,1);
func=@(s) 1/h+(pi^2*h*s.^2)/4;
b(N)=quad8(func,0,1);
%c=a;
a=[0,a];
c=[a,0];


%得到总载荷向量Total load vector；
function d=TLV(x,d,h,N)
format long
for i=1:N-1
    %被积函数为f(x(i)+h*s)*h*s + f(x(i+1)+h*s)*h*(1-s)
    %func=@(s) pi^2*(sin(pi*(x(i)+h*s)./2)+sin(pi*(x(i+1)+h*s)./2))/2;
    %func=@(s) pi^2*h*(s.*sin((pi*(x(i)+h*s))/2)+(1-s).*sin((pi*(x(i+1)+h*s))/2)/2;
    func=@(s) pi^2*h*s.*sin((pi*(x(i)+s.*h))/2)/2+pi^2*h*s.*sin((pi*(x(i+1)+s.*h))/2)/2;
    d(i)=quad(func,0,1);
end
func=@(s) pi^2*h*s.*sin((pi*(x(N)+s.*h))/2)/2;
d(N)=quad(func,0,1);


%计算有限元的主程序；
function [u,v]=compute_ODFE(minx,maxx,N)

%计算步长；
h=(maxx-minx)/N;

%初始化a,b,c,d,u,v;
a=zeros(1,N-1);
b=zeros(1,N);
c=zeros(1,N-1);
d=zeros(1,N);
u=zeros(1,N);
v=zeros(1,N+1);

%计算小区间上的单元点对应的x值；
for i=1:N+1
    x(i)=minx+(i-1)*h;
    %顺便求得精确解；
    v(i)=sin(pi*x(i)/2);
end

%先得到总刚度矩阵的对角线上的元素；
[a,b,c]=TSM(x,a,b,c,h,N);
d=TLV(x,d,h,N);

%现在再解方程组得到数值解；
u=threedia(a,b,c,d');
%直接用右除来求解方程组；
%K=diag(a,-1)+diag(b)+diag(c,1);
%u=Kd';
%u=[0,u'];
u=[0,u];

%计算误差；
e=abs(u-v);

%绘制图形；
figure;
subplot(1,3,1);
plot(x,u,'.');
title('数值解');
subplot(1,3,2);
plot(x,v,'r');
title('精确解');
subplot(1,3,3);
plot(x,e,'-');
title('误差');


%追赶法求解线性方程组；
function x=threedia(a,b,c,f)
%示例系数矩阵：【2，-1，0，0    右边项：6
%                -1，3，-2，0           1
%                0，-1，2，-1           0
%                0，0，-3，5】          1
%参数按照这种方式输入：a=[0,-1,-1,-3];b=[2,3,2,5];c=[-1,-2,-1,0];f=[6,1,0,1]';
N=length(f);
x=zeros(1,N);
y=zeros(1,N);
beta=zeros(1,N);
gramma=zeros(1,N);
beta(1)=b(1);
for i=1:N-1
    gramma(i)=c(i)/beta(i);
    beta(i+1)=b(i+1)-a(i+1)*gramma(i);
end
%追的过程
y(1)=f(1)/beta(1);
for i=2:N
    y(i)=(f(i)-a(i)*y(i-1))/beta(i);
end
%赶的过程
x(N)=y(N);
for i=N-1:-1:1
    x(i)=y(i)-gramma(i)*x(i+1);
end

-------------------------------------------补充更新-------------------------------------------

这个本来已经发表出去了的，但是我回过头去看的时候写的太匆忙了，现在想要改一改整理好重新上传一次。然后还有一位同学私信我说最近两天没有更新了，很抱歉，这两天可能是学习的劲头有点不足，又加之有点其他的事情所以一直没有学了，然后想修改这一篇的时候费了好大劲重新理清思路，学习还是要持续啊，不然浪费那么多时间回顾之前的也是太心累了，一片混乱的感觉。

我主要改了一下的就是追赶法的数值算法，是直接复制了CSDN博客上某为博友的追赶法的代码，（有现成的当然是直接拿来用啦~也怪自己以前本科的时候不懂学习方法，那时候数值分析学的现在都是基础啊，都可有用了！）另外，关于那个自动调用函数的我不打算改啦，下次如果遇到了就都先自己手算一下再修改代码里面函数部分就好啦，我把方程组的元素也放在代码里面注释好了，只要按照注释代入一下自变量即可，需要注意的地方就是编的时候s有几个给地方是.*的形式，以下是追赶法的结果，看起来和直接用右除的结果是一样的，自己看需要决定啦，也不是所有情况两个算出来的结果都是一样的。个人认为第二份代码比较整齐。

代码如下：

%追赶法求解三对角线性方程组，Ax=b，A用一维数组a，c，d存储。
function x=crout(a,c,d,b)
%数组a存储三角矩阵A的主对角线元素，c、d存储主对角线上边下边带宽为1的元素
n=length(a);
n1=length(c);
n2=length(d);
%错误检查
if n1~=n2%存储矩阵的数组维数错误
    error('MATLAB:Crout:不是三对角矩阵，参数数组中元素个数错误.');
elseif n~=n1+1
    error('MATLAB:Crout:不是三对角矩阵，参数数组中元素个数错误.');
end
%初始化
L=zeros(n);%生成n*n的全零矩阵
U=zeros(n);
p=1:n;
q=1:n-1;
x=1:n;
y=1:n;
%追赶法程序主体
p(1)=a(1);
for i=1:n-1
    q(i)=c(i)/p(i);
    p(i+1)=a(i+1)-d(i)*q(i);%d的下标改为1到n-1
end
%正解y
y(1)=b(1)/p(1);%用x存储y
for i=2:n
    y(i)=(b(i)-d(i-1)*y(i-1))/p(i);
end
%倒解x
x(n)=y(n);
for i=(n-1):-1:1
    x(i)=y(i)-q(i)*x(i+1);
end

%计算有限元的主程序；
function [u,v]=compute_ODFE(minx,maxx,N)

%计算步长；
h=(maxx-minx)/N;

%初始化a,b,c,d,u,v;
a=zeros(1,N-1);
b=zeros(1,N);
c=zeros(1,N-1);
d=zeros(1,N);
u=zeros(1,N);
v=zeros(1,N+1);

%计算小区间上的单元点对应的x值；
for i=1:N+1
    x(i)=minx+(i-1)*h;
    %顺便求得精确解；
    v(i)=sin(pi*x(i)/2);
end

%先得到总刚度矩阵的对角线上的元素；
[a,b,c]=TSM(x,a,b,c,h,N);
%a为-1对角线上的元素；
%b为主对角线上的元素；
%c为1对角线上的元素；
d=TLV(x,d,h,N);
%方程组的右端项；

%现在再解方程组得到数值解；
u=crout(b,c,a,d);
%直接用右除来求解方程组；
%K=diag(a,-1)+diag(b)+diag(c,1);
%u=Kd';
u=[0,u];

%计算误差；
e=abs(u-v);

%绘制图形；
figure;
subplot(1,3,1);
plot(x,u,'.');
title('数值解');
subplot(1,3,2);
plot(x,v,'r');
title('精确解');
subplot(1,3,3);
plot(x,e,'-');
title('误差');

%得到总载荷向量Total load vector；
function d=TLV(x,d,h,N)
format long
for i=1:N-1
    %被积函数为f(x(i)+h*s)*h*s + f(x(i+1)+h*s)*h*(1-s)
    %直接计算之后的结果；
    func=@(s) pi^2*h*s.*sin((pi*(x(i)+s.*h))/2)/2+pi^2*h*s.*sin((pi*(x(i+1)+s.*h))/2)/2;
    d(i)=quad(func,0,1);
end
%最后一项被积分函数表示：f(x(N)+h*s)*h*s
%直接计算之后的结果；
func=@(s) pi^2*h*s.*sin((pi*(x(N)+s.*h))/2)/2;
d(N)=quad(func,0,1);


%得到总刚度矩阵Total stiffness matrix的三条对角线上的元素；
%a为-1对角线上的元素；
%b为主对角线上的元素；
%c为1对角线上的元素；
%h为步长；
function [a,b,c]=TSM(x,a,b,c,h,N);
format long
%计算对角线上的元素值；
for i=1:N-1
    %先计算主对角线上的元素；
    %被积函数表示：p(x(i)+h*s)/h+q(x(i)+h*s)*h*s^2 + p(x(i+1)+h*s)/h+q(x(i+1)+h*s)*h*(1-s)^2 ;
    %还需要额外再算最后一项的积分；
   
    %再计算-1对角线上的元素：
    %被积分函数表示：q(x(i+1)+h*s)*h*s*(1-s)-p(x(i+1)+h*s)/h;
  
    %这里是直接计算好函数之后直接的结果主对角线上的元素；
    func=@(s) 2/h+((pi^2)*h*(s.^2+(1-s).^2))/4;
    b(i)=quad8(func,0,1);
    %这里是直接计算好函数之后直接的结果-1对角线上的元素；
    func=@(s) (pi^2*h*s.*(1-s))/4-1/h;
    a(i)=quad8(func,0,1);
end

%最后一项被积分函数表示：p(x(N)+h*s)/h+q(x(N)+h*s)*h*s^2;
%这里是直接计算好函数之后直接的结果；
func=@(s) 1/h+(pi^2*h*s.^2)/4;
b(N)=quad8(func,0,1);
c=a;