边值问题的研究.docx

西北农林科技大学理学院应用数学系微分方程数值解结课论文论文题目 边值问题的研究2016年 1 月 14 日一问题重述对于下列边值问题：其中A为学号的倒数第2位，B为学号的倒数第1位。(1)差分：截断误差、稳定性、收敛半径、递推（隐式）或方程组（显式）(2)有限元：刚度矩阵、算法步骤及代码二问题分析题目明确指出使用差分方法和有限元解法。什么都不管先构造一种差分格式，然后对求解区域做划分将问题离散化，从微分方程的定解问题转化为求线性代数方程的解，以便于能够使用计算机进行计算。在这里选用的是中心差分法，同时将边界进行处理，同时用Ritz有限元法和Galerkin法有限元法尝试去得到结果，最后再去比较两种解法所得到结果的精确性，分析相容性和截断误差等等。三解题过程1首先建立差分格式，考虑两点的边值问题，由题目知道建立中心差分格式如下对求解区间做网格划分，在a到b之间取N+1个节点，定义为xi（i取1到N）即将区间I=a,b分为N个小区间由此得到区间的一个网格剖分。记。用表示网格内点，的集合，表示内点和界点的集合。取相邻节点的中点，称为半整数点。由节点又构成的一个对偶剖分。用差商代替微商，将方程（1.1）在内点离散化. 逼近边值问题（1.1）（1.2）的差分方程为：当网格均匀，即时差分方程简化为这相当于用一阶中心差商，二阶中心差商依次代替（1.1）的一阶微商和二阶微商的结果。这个方程就是中心差分格式。式（1.4）用方程组展开：这是一个以为未知量的线性方程组。到此为止，中心差分格式展开完毕，接下来处理方程(1.1)将方程在节点离散化，由泰勒公式展开得：所以截断误差为下一步是分析差分格式的稳定性差分格式的截断误差：，而边界条件的截断误差为收敛性和稳定性是从不同角度讨论差分法的精确情况，稳定性主要是讨论初值的误差和计算中的舍入误差对计算结果的影响，收敛性则主要讨论推算公式引入的截断误差对计算结果的影响.使用既收敛有稳定的差分格式才有比较可靠的计算结果，这也是讨论收敛性和稳定性的重要意义.截断误差：，即。差分方程组的解满足：其中a、b代表边界点，代表边界点的取值。上式给出了差分方程的解的误差估计，而且表明当差分解收敛到原边值问题的解，收敛速度为。2接下来是有限元的解法从Ritz法出发，单元刚度矩阵为：按规则组装成总刚度矩阵。令其中以及则有限元方程为从Galerkin有限元法出发，Galerkin有限元方程为：系数矩阵第j行只有三个非零元素，即这里第一行只有两个非零元素：第n行只有两个非零元素：和方程的右端项四求解过程其精确解为。算例中。(1)从Ritz法出发以将积分区间等分为10份为例，则步长，记为。 为：以步长取为h=1/10为例，从Ritz法出发的有限元法得到的数值解与精确解为Ritz数值解精确解001.15401.11352.22452.14803.20253.09454.07903.94404.84504.68755.49155.31606.00955.82056.39006.19206.62406.42156.702506.5000图像为分析：最大误差为0.202500，Ritz有限元法求解两点边值问题很接近精确解。以步长取为h=1/50为例，从Ritz法出发的有限元法得到的数值解与精确解图像为最大误差为0.002025步长1/101/501/1001/500Ritz最大误差0.20250.0441000.020250.004491分析：最大误差为0.02025，Ritz有限元法求解两点边值问题很接近精确解，且步长越大，误差越小。（2）从Galerkin法出发以将积分区间等分为10份为例，则步长，记为。 为：Galerkin有限元法最大误差:0.815000，图像为： 以将积分区间等分为100份为例，图像为：分析：最大误差为0.080150，Galerkin有限元法求解两点边值问题很接近精确解，且步长越大，误差越小。步长1/101/501/1001/500Calerkin最大误差0.801500.160060.0801500.016006最后收敛性和误差分析：令和分别表示精确解和有限元解在剖分区间节点处的值，收敛性表示为记最大误差为err,则问题转化为在方程中，已知h和err，求解M和k的拟合问题。在Matlab中拟合采用最小二乘法实现。对和进行最小二乘幂函数拟合，求得从Ritz法出发的误差阶为k=0.9612,M=0.4115.对和进行最小二乘幂函数拟合，求得从Galerkin法出发的误差阶为k=1.004,M=3.061.五操作代码主程序：function Ritz(a,b,N)% -D2u=a*x+b，u(0)=0，Du(1)=0%a=9;b=7; %a为学号倒数第二位，b为学号倒数第一位，%N为剖分份分数% 调用格式：Ritz(9,7,10); %N为剖分份分数syms x;ua=0; %区间界点ub=1; %区间界点u_a=0; %左边界条件du_b=0; %右边界条件p=1;q=0;f=a*x+b; %f=9x+7h=1/N;X=0:h:1; K=Stiff_matrix(p,q,h,N,X); %得到总刚度矩阵KB=vector(p,q,X,h,N,f); %得到BU = KB; %差分解 %处理右边值条件u_b = (2*h*du_b-U(end-1)+4*U(end)/3;U=0;U;u_b; %精确解y0 = dsolve(-D2y=9*X+7,y(0)=0,Dy(1)=0,X);precise_value=double(subs(y0); deta=U-precise_value; deta_max=max(abs(deta); %最大误差fprintf(最大的误差是%fn,deta_max)% 差分解与精确解对比表figure;subplot(1,2,1);plot(X,U,b*,X,precise_value,r-);xlabel(x);ylabel(u);title(差分解与精确解对比图);legend(差分解,精确解,Location,best);% 误差图subplot(1,2,2);plot(X,deta);xlabel(x);ylabel(u);end子程序：得到刚度矩阵K：function K=Stiff_matrix(p,q,h,N,X) syms x;K=zeros(N-1,N-1); diag_0=zeros(N-1,1); %确定K的对角元diag_1=zeros(N-2,1); %确定K的偏离对角的上对角元diag_2=zeros(N-2,1); %确定K的偏离对角的下对角元 % 获取对角元素for j=1:N-1 diag_0(j)=(subs(p,x,(X(j)+X(j+1)/2)+(subs(p,x,(X(j)+X(j+1)/2+h)+(subs(q,x,X(j+1)*(h2);end% 获取A的第三条对角for j=1:N-2 diag_2(j)=-(subs(p,x,(X(j+1)+X(j+2)/2)-(subs(0,x,X(j+2)*h)/2;end%获取A的第二条对角for j=1:N-2 diag_1(j)=-(subs(p,x,(X(j+1)+X(j+2)/2)+(subs(0,x,X(j+1)*h)/2;endK=diag(diag_0)+diag(diag_1,1)+diag(diag_2,-1); K(N-1,N-1)=2;K(N-1,N-2)=-2;得到B：function B=vector(p,q,x,h,N,f)B=zeros(N-1,1); syms x;X=0:h:1;for j=1:N-1; B(j)=(h2)*(subs(f,x,X(j+1);end%系数矩阵B(1)=B(1)+0*(subs(p,x,(X(2)+X(1)/2)+(subs(0,x,X(2)*h/2);B(N-1)=3*B(N-1)+2*0*(subs(p,x,(X(N)+X(N+1)/2)-(subs(0,x,X(N)*h/2)主程序：p=1;q=0;a = 9;b = 7;N = 100;%剖分份数h=1/N;x=0:h:1;A_Galerkin=matirix(a,b,p,q,h,N);values_f_Galerkin=vector1(a,b,x,h,N);U_Galerkin=A_Galerkinvalues_f_Galerkin; y0 = dsolve(-D2y=a*x+b,y(0)=0,Dy(1)=0,x);precise_value=double(subs(y0);% Galerkin有限元法解与精确解对比图figure;%subplot(1,2,1);plot(x,0;U_Galerkin,b.-,x,precise_value,r.-);xlabel(x);ylabel(u);title(Galerkin有限元法解与精确解对比图);legend(Galerkin数值解,精确解,Location,best);err_Galerkin=max(abs(0;U_Galerkin-precise_value);fprintf(sprintf(Galerkin有限元法最大误差:%fn,err_Galerkin);子程序：% Galerkin有限元法求解一维问题%构造系数矩阵function A_Galerkin=matirix(a,b,p,q,h,N),A_Galerkin=zeros(N);for i=3:N fun1_Galerkin=(ks)-p./h+h.*q.*ks.*(1-ks); fun2_Galerkin=(ks)p./h+h.*q.*(ks.2)+p./h+h.*q.*(1-ks).2); fun3_Galerkin=(ks)-p./h+h.*q.*ks.*(1-ks); A_Galerkin(i-1,i-2)=integral(fun1_Galerkin,0,1); A_Galerkin(i-1,i-1)=integral(fun2_Galerkin,0,1); A_Galerkin(i-1,i) =integral(fun3_Galerkin,0,1);endfun2_Galerkin=(ks)p./h+h.*q.*(ks.2)+p./h+h.*q.*(1-ks).2);A_Galerkin(1,1)=integral(fun2_Galerkin,0,1);fun3_Galerkin=(ks)-p./h+h.*q.*ks.*(1-ks);A_Galerkin(1,2)=integral(fun3_Galerkin,0,1);fun3_Galerkin=(ks)-p./h+h.*q.*ks.*(1-ks);A_Galerkin(N,N-1)=integral(fun3_Galerkin,0,1);fun4_Galerki