您的当前位置:首页正文

两点边值问题的差分求解

2022-12-29 来源:九壹网
微分方程数值解实验报告

姓名 丁建伟 学号 200708020211 日期 2010.11.25 徐强 实验项目 两点边值问题的差分求解 指导教师 一、上机实验的问题和要求(需求分析): 实验内容: (I) 分别在步长h=1/20,1/40,1/80,1/160情形下用中心差分格式计算齐次两点边值问题-u"=f,u(0)=u(1)=0。其中f(x) = 100*exp(-10*x),精确解为u(x) = 1 - (1-exp(-10))*x - exp(-10*x) (II) 给出差分解近似精确解在无穷范数和L2范数下的误差阶。 目的与要求: 掌握中心差分格式的程序实现 掌握分析算法误差的方法 二、程序设计的基本思想,原理和算法描述: 基本思想及原理: 做均匀网格剖分: 0x0x1xN1,分点xiih在节点1h步长n xd2uf(x) i处,对微分方程离散化2dxu(xi1)2u(xi)u(xi1)d2uh2d4u3 O(h)224hdxi12dxiu(xi1)2u(xi)u(xi1) f(xi)Ri(u) 有2hh2d4u3R(u)O(h) 4 其中 i12dxi记u在节点xk,k0~N数值解为 uk,k0~N, ui12uiui1fi, 2h则有Lhui:比较知 所以Lhu(xi):f(xi)Ri(u) Ri(u)Lhu(xi)Lui Lh代替微分算子L产生的误差 表示用差分算子2O(h)。 称之为(局部)截断误差。这里关于h的阶为 注意所以Luif(xi) Ri(u)Lhu(xi)f(xi) 由此知:(局部)截断误差可视为差分格式,将数值解换成相应真解值后,左端减右端,再做Taylor展式获得的(可作为计算公式)。 方程的联立形式(中心差分格式) ui12uiui1fi,i1~n12h u0,u00N矩阵形式 AUb(其中 A 是三对角矩阵) 又因为A是三对称矩阵,而且符合追赶法的使用条件,故可用追赶法求解U的解。 三、主要程序代码或命令: #include #include #define MAX 200 /*预定义数组大小*/ void main() { int n,i; /*初始化阶数n*/ float u[MAX],y[MAX]; float F[MAX],f[MAX],m[MAX]; float h,x; /*步长和剖分点*/ printf(\"请输入等分数n值:\"); scanf(\"%d\ /*读入阶数*/ h=1/float(n); m[1]=-0.5; /*使用追赶法求解系数矩阵三对称的线性方程组*/ for(i=2;i<=n-2;i++) m[i]=-1/(2+m[i-1]); for(i=1,x=h;x<1.0;x+=h,i++) {F[i]=100*exp(-10*x); f[i]=1-(1-exp(-10))*x-exp(-10*x); } y[1]=F[1]/(2/(h*h)); for(i=2;i<=n-1;i++) y[i]=(F[i]+(1/(h*h)*y[i-1]))/(2/(h*h)+(1/(h*h)*m[i-1])); u[n-1]=y[n-1]; for(i=n-2;i>=1;i--) u[i]=y[i]-m[i]*u[i+1]; for(i=1;i

因篇幅问题不能全部显示,请点此查看更多更全内容