中心差分解两点边值问题_第1页
中心差分解两点边值问题_第2页
中心差分解两点边值问题_第3页
中心差分解两点边值问题_第4页
中心差分解两点边值问题_第5页
已阅读5页,还剩5页未读 继续免费阅读

下载本文档

版权说明:本文档由用户提供并上传,收益归属内容提供方,若内容存在侵权,请进行举报或认领

文档简介

1、一题目用中心差分格式计算如下两点边值问题已知其精确解为二理论作为模型,考虑两点边值问题:1.11.2假定是给定的常数。1. 建立差分格式1.区域网格剖分首先取个节点:将区间分成个小区间:于是得到区间的一个网格剖分。记,称为网格最大步长。用表示网格内点,的集合,表示内点和界点的集合。取相邻节点的中点,称为半整数点。则由节点又构成的一个网格剖分,称为对偶剖分。2.微分方程的离散,建立相应差分格式用差商代替微商,将方程1.1在内点离散化.注意对充分光滑的,由Taylor展式有1.3 1.5由1.5减1.4,并除以,得1.6令则由1.31.6知,边值问题的解满足方程: 1.7 其中 1.8为差分算子的

2、截断误差,舍去,便得逼近边值问题1.11.2的差分方程: 1.9i=1,2,N-1,由方程1.71.9,截断误差可表示为 1.10当网格均匀,即时差分方程1.9简化为 1.11这相当于用一阶中心差商,二阶中心差商依次代替1.1的一阶微商和二阶微商的结果。这个方程就是中心差分格式。截断误差为: 1.12所以截断误差按或的阶为。在此题中, ,因为r=0方程1.11的系数对角矩阵是三对角矩阵。我们可以用消元法或迭代法求解方程组1.11.2式1.11用方程组展开:写成矩阵形式为:2.收敛性分析根据1.10我们引进误差则误差函数满足以下差分方程:于是收敛性及收敛速度的估计问题,就归结到通过右端截断误差估

3、计误差函数的问题。由1.12我们知,有从而差分方程满足相容条件。假设引进记号,设则可将1.9改写为将差分解表成 2.1其中满足 2.2而满足 2.3先估计,由 2.4据差分格林公式再利用柯西不等式,有常数使 2.5将不等式2.6用于2.5右端,则 2.6解差分方程2.2,易得从而这样, 2.7利用范数,从2.7推出 2.8因为因此 2.9联结2.12.7及2.9即得差分解的先验估计: 2.10其中不等式2.10说明差分解连续依赖于右端和边值,因此差分格式1.11关于右端及边值稳定.根据定理1.1 : 假设边值问题的解u充分光滑,差分方程按满足相容条件且关于右端稳定,则差分解按收敛到边值问题的解

4、,且有和相同的收敛阶。所以差分方程的解的收敛速度为。三程序代码:clcclfclfsyms x;a=1; %区间界点b=2; %区间界点p=exp(x); %这是p函数q=sin(x)+1+x; %这是q函数f=-exp(x)*(2*x+1)+(sin(x)+1+x)*x*(x-1);%这是f函数r=0; %这是r函数.N=10; %将区间划分的等分,这里控制!h=(b-a)/N; %这里确定步长value_of_f=zeros(N-1,1);%这是fdiag_0=zeros(N-1,1);%确定A的对角元diag_1=zeros(N-2,1);%确定A的偏离对角的上对角元diag_2=zer

5、os(N-2,1);%确定A的偏离对角的下对角元X=a:h:b;u_a=0; %边界条件u_b=2; %边界条件for j=2:N diag_0(j-1)=(subs(p,x,(X(j+1)+X(j)/2)+(subs(p,x,(X(j-1)+X(j)/2)/(h2)+(subs(q,x,X(j);end %获取对角元素for j=3:N diag_2(j-2)=-(subs(p,x,(X(j-1)+X(j)/2)/(h2)-subs(r,x,X(j)/(2*h);end %获取A的第三条对角for j=2:N-1 diag_1(j-1)=-(subs(p,x,(X(j+1)+X(j)/2)/

6、(h2)+subs(r,x,X(j)/(2*h);end %获取A的第二条对角for j=2:N; value_of_f(j-1)=subs(f,x,X(j);end %获取F值value_of_f(1)=value_of_f(1)+u_a*(subs(p,x,(X(2)+X(1)/2)/(h2);value_of_f(N-1)=value_of_f(N-1)+u_b*(subs(p,x,(X(N)+X(N+1)/2)/(h2);A=diag(diag_0)+diag(diag_1,1)+diag(diag_2,-1);%组装系数矩阵format longU=inv(A)*value_of_f

7、 %差分解%fprintf('%11.5f',U)fprintf('n');dx=X(2:N);precise_value=dx.*(dx-1) %精确解%fprintf('%11.5f',precise_value)deta=U-precise_value' ; %误差deta_max=max(abs(deta);%最大误差fprintf('最大的误差是%fn',deta_max)plot(X(2:N),U,'b*',X(2:N),precise_value,'r-') %差分解与精确解比

8、照表figure();plot(X(2:N),deta) %误差图结果: X的值步长h2.12.22.32.42.52.62.72.82.9最大误差0.10.110150.240260.390330.560370.750380.960381.190311.440221.710120.0003800.050.110030.240060.390080.560090.750090.960081.190071.440051.710030.0000950.0250.110000.240010.390020.560020.750020.960021.190021.440011.710010.0000240.01250.110000.240000.390000.560010.750010.960001.190001.440001.710000.000006精确解0.110000.240000.390000.560000.750000.960001.190001.440001.710000以下仅给出步长为N=20,h=0.05的精确值和差分值图与误差图,其他的只要修改程序中的步

温馨提示

  • 1. 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。图纸软件为CAD,CAXA,PROE,UG,SolidWorks等.压缩文件请下载最新的WinRAR软件解压。
  • 2. 本站的文档不包含任何第三方提供的附件图纸等,如果需要附件,请联系上传者。文件的所有权益归上传用户所有。
  • 3. 本站RAR压缩包中若带图纸,网页内容里面会有图纸预览,若没有图纸预览就没有图纸。
  • 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
  • 5. 人人文库网仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对用户上传分享的文档内容本身不做任何修改或编辑,并不能对任何下载内容负责。
  • 6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
  • 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。

评论

0/150

提交评论