微分方程数值解实习题

数学课程微分方程求解练习题及答案微分方程是数学中非常重要的一门课程，它在许多科学领域中有着广泛的应用。

为了更好地掌握微分方程的解题技巧，下面将给出一些微分方程求解的练习题及其答案。

练习一：一阶线性微分方程1. 求解微分方程：dy/dx + y = 2x解答：首先将该微分方程转化为标准形式：dy/dx = 2x - y然后可以使用分离变量的方法进行求解，将变量分离得到：dy/(2x - y) = dx对等式两边同时积分，得到：∫(1/(2x - y))dy = ∫dx通过对右边的积分，得到：ln|2x - y| = x + C1 (其中C1是常数)将等式两边取e的指数，得到：2x - y = Ce^x其中C = e^C1是一个任意常数，所以方程的通解为：y = 2x - Ce^x （其中C为常数）2. 求解微分方程：dy/dx + 2y = e^x解答：将该微分方程转化为标准形式：dy/dx = e^x - 2y然后使用分离变量的方法进行求解，得到：dy/(e^x - 2y) = dx对等式两边同时积分，得到：∫(1/(e^x - 2y))dy = ∫dx通过对右边的积分，得到：(1/2)ln|e^x - 2y| = x + C2 （其中C2是常数）再次将等式两边取e的指数，得到：e^x - 2y = Ce^2x其中C = e^C2是一个任意常数，所以方程的通解为：y = (1/2)e^x - (C/2)e^2x （其中C为常数）练习二：二阶微分方程1. 求解微分方程：d^2y/dx^2 + 4dy/dx + 4y = 0解答：首先将该微分方程的特征方程写出来：r^2 + 4r + 4 = 0解特征方程，得到特征根为：r = -2由于特征根为重根，所以方程的通解形式为：y = (C1 + C2x)e^(-2x) （其中C1和C2为常数）2. 求解微分方程：d^2y/dx^2 + dy/dx - 2y = 0解答：首先将该微分方程的特征方程写出来：r^2 + r - 2 = 0解特征方程，得到特征根为：r1 = 1，r2 = -2所以方程的通解形式为：y = C1e^x + C2e^(-2x) （其中C1和C2为常数）这里给出了一些微分方程求解的练习题及其答案，通过练习这些题目，相信可以增强对微分方程的理解和掌握。

微分方程练习练习1 求解范德堡（vander pol）方程练习2 单摆运动图4.3中一根长的细线，一端固定，另一端悬挂质量为的小球，在重力作用下，小球处于竖直的平衡位置. 现使小球偏离平衡位置一个小的角度，然后使其自由运动，在不考虑空气阻力情形下，小球将沿弧线作周期一定的简谐运动.为平衡位置，在小球摆动过程中，当与平衡位置夹角为时，小球所受重力在其动运轨迹的分量为（负号表示力的方向使减少），由牛顿第二定律可得微分方程（4.12）设小球初始偏离角度为，且初速为0，式（4.12）的初始条件为（4.13）当不大时，，式（4.12）化为线性常系数微分方程图4.3（4.14）解得（4.15）简谐运动的周期为.现在的问题是：当较大时，仍用近似，误差太大，式（4.12）又无解析解，试用数值方法在两种情况下求解，画出的图形，与近似解（4.15）比较，这里设.练习3 捕食与被捕食当鲨鱼捕食小鱼，简记为乙捕食甲，在时刻，小鱼的数量为，鲨鱼的数量为，当甲独立生存时它的（相对）增长率与种群数量成正比，即有，为增长率，而乙的存在使甲的增长率减少，设减少率与乙的数量成正比，而得微分方程（4.16）比例系数反映捕食者掠取食饵的能力.乙离开甲无法生存，设乙独自存在时死亡率为，，甲为乙提供食物，使乙的死亡率降低，而促其数量增长，这一作用与甲的数量成正比，于是满足（4.17）比例系数反映甲对乙的供养能力，设若甲，乙的初始数量分别为（4.18）则微分方程（4.16），（4.17）及初始条件（4.18）确定了甲，乙数量、随时间变化而演变的过程，但该方程无解析解，试用数值解讨论以下问题：（1）设，求方程（4.16），（4.17）在条件（4.18）下的数值解，画出的图形及相图，观察解的周期变化，近似确定解的周期和的最大、小值，近似计算在一个周期内的平均值.（2）从式（4.16）和（4.17）消去得到（4.19）解方程（4.19），得到的解即为相轨线，说明这是封闭曲线，即解确为周期函数.（3）将方程（4.17）改写为（4.20）在一个周期内积分，得到一周期内的平均值，类似可得一周期内的平均值，将近似计算的结果与理论值比较.进一步练习（1）编写改进欧拉公式求微分方程数值解的程序，并用其与ode23求下列微分方程数值解，对二者作出比较.a）或.b)（Bessel方程，这里令，其精确解为）.c).（2）倒圆锥形容器，上底面直径为1.2m，容器的高亦为1.2m，在锥尖的地方开有一直径为3cm的小孔，容器装满水后，下方小孔开启，由水利学知识可知当水面高度为时，水从小孔中流出的速度为为重力加速度，若孔口收缩系数为0.6（即若一个面积单位的小孔向外出水时，水柱截面积为0.6），问水从小孔中流完需多少时间？2分钟时，水面高度是多少？（3）一只小船渡过宽为的河流，目标是起点正对着的另一岸上点，已知河水流速与船在静水中的速度之比为.（a）建立小船航线的方程，求其解析解.（b）设，用数值解法求渡河所需时间，任意时刻小船的位置及航行曲线，作图并与解析解比较.（c）若流速为0,0.5,2 (m/s)，结果将如何？（4）研究种群竞争模型. 当甲、乙两个种群各自生存时，数量演变服从下面规律其中，分别为时刻甲，乙两个种群的数量，为其固有增长率，为它们的最大容量，而当这两个种群在同一环境中生存时，由于乙消耗有限资源而对甲的增长产生影响，将甲的方程修改为（4.22）这里的含意是：对于供养甲的资源而言，单位数量乙（相对）的消耗率为单位数量甲（相对）消耗的倍，类似地，甲的存在亦影响乙的增长，乙的方程应改为（4.23）给定种群的初始值为（4.24）及参数后，方程（4.22）与（4.23）确定了两种群的变化规律，因其解析解不存在，试用数值解法研究以下问题：（a）设，计算，画出它们的图形及相图，说明时间充分大以后的变化趋势（人们今天看到的已经是自然界长期演变的结果）.（b）改变，但与不变（保持），分析所得结果，若，再分析结果，由此你得到什么结论，请用各参数生态学上的含义作出解释.（c）试验当时会有什么结果；当时又会出现什么结果，能解释这些结果吗？。

常微分方程习题 《李立康》习题1.用Euler 方法求初值问题⎩⎨⎧=-='0)0(21u tuu 在1=t 时的近似解（取41=h ）。

2.初值问题1300u u u()⎧⎪'=⎨⎪=⎩ 有解3223/u(t )t ⎛⎫= ⎪⎝⎭。

但若用Euler 方法求解，对一切N T ,和HTh =,都只能得到N t u t , (2)1,0==，试解释此现象产生的原因。

3.用Euler 方法计算⎩⎨⎧=='1)0(u uu 在1=t 处的值，取161和41=h ，将计算结果与精确值e =）1（u 相比较。

4.设),(u t f 满足定理2.1的条件，对改进Euler 法（2.10）式证明： （1）其局部截断误差为)()(1243h O t u h -'''-；（2）当1<hL 时，其整体截断误差满足：)1(22--≤Lt n lT m e hLRe εε （3）方法具有二阶收敛速度且稳定。

5.导出用改进Euler 法求解⎩⎨⎧=='1)0(u uu 计算公式mmh h u ⎪⎪⎭⎫ ⎝⎛-+=22 取41=h 计算)1(u 的近似值，并与习题3的结果比较。

6.就初值问题⎩⎨⎧=+='0)0(u bat u 分别导出用Euler 方法和改进Euler 法求近似解的表达式，并与真解bt t au +=22相比较。

7.证明改进Euler 法的绝对稳定区域是整个左半平面0)Re(<h 。

8.对初值问题⎩⎨⎧=-='1)0(2u u u 用41=h 的Euler 方法求解，求出实际计算值t u 与真解tu +=11在)1(u 处的误差，并将它与定理2.3的估计式（2.22）式相比较。

9.证明：Runge-Kutta 方法中);,(h u t ϕ关于u 或t 满足Lipschitz 条件的充分条件是),(u t f 关于t 或u 满足Lipschitz 条件。

《微分方程数值解法》复习、练习题第一章复习题1、建立差分格式的三个主要步骤（三个离散化）。

2、差分格式的相容性、收敛性概念。

3、Poisson 方程的5点菱形差分格式，矩形、非矩形区域情形边界条件的处理（离散化）。

4、对长方形区域作正方形网格剖分，求解Poisson 方程边值问题的五点菱形差分格式，按什么顺序对节点编号，可使差分方程带宽更窄？（按短方向排）5、差分方程有哪些共同特性，求解选用哪类方法？（大型稀疏，带状，主对角占优等，一般采用迭代法）多重网格等略。

6、极值原理。

7、5点菱形差分格式求解Poisson 方程第一边值问题的收敛性。

第一章练习题1、设有边值问题=?+??-=-==<<<<=?====x u n u u y u u y x x u y y x x 2,1122.00,3.00,2.003.00取h =0.1的正方形网格。

（1）用5点菱形格式在内点建立差分格式；（2）用截断误差为)(2h O 的方法离散化第三边界条件（有两种方式）；（3）写出整理后的差分方程的矩阵形式=??????? ????????? ?D C B A u u u u2、定义方形算子如下：(),1,11,11,11,1,2142i j i j i j i j i j i j u u u u u u h---++-++=+++- 试讨论5点方形差分方程,,i j i j u f =逼近微分方程(,)u f x y ?=的截断误差是几阶？3、设有{}220,(,)0,1ln (1)u x y x y u x y ?Ω?=∈Ω=<，取h =1/3，列出5点方形差分格式所得的差分方程。

第二章复习题1、差分格式稳定性与收敛性的定义。

2、有关求特征值的几个结论。

3、判断稳定性的矩阵法和Fourier 分析法（Von-Neumann 条件）的应用。

4、显隐格式在一般情况下的优缺点。

5、熟悉古典显、隐格式，六点对称隐格式（C-N 格式）。

微分⽅程数值解法答案包括基本概念，差分格式的构造、截断误差和稳定性，这些内容是贯穿整个教材的主线。

解答问题关键在过程，能够显⽰出你已经掌握了书上的内容，知道了解题⽅法。

这次考试题⽬的类型：20分的选择题，主要是基本概念的理解，后⾯有五个⼤题，包括差分格式的构造、截断误差和稳定性。

习题⼀1．略2． y y x f -=),(，梯形公式：n n n n n n y hh y y y h y y )121(),(2111+-+=+-=+++，所以0122)1(01])121[()121()121(y hh y h h y h h y hhn h h n n n +--+--+-+=+-+==+-+= ，当0→h 时，x n e y -→。

同理可以证明预报-校正法收敛到微分⽅程的解.3．局部截断误差的推导同欧拉公式;整体截断误差:++++++-++≤1),())(,(11111n nx x n n n n n n n dx y x f x y x f R εε11)(++-++≤n n n y x y Lh R ε，这⾥R R n ≤ ⽽111)(+++-=n n n y x y ε，所以 R Lh n n +=-+εε1)1(，不妨设1()]11111[1111101---++-+-+-≤≤-+-=n n n n Lh Lh Lh R Lh Lh R Lh εεε ]1[2)(02)(00-+≤--x X L x X L eLh R eε4．中点公式的局部截断误差： dx x y x f hx y h x f x y x f yx y n n x x n n n n n n))](,(2)(,2())(,([)(11*1?+++-=-++dx x y x f hx y h x f h x y h x f h x y x y dxx y x f hx y h x f hx y h x f h x y h x f x y x f n n n n x x n n n n n n n x x n n n n n n n n))](,(2)(,2())2(,2([)]2()([))](,(2)(,2())2(,2())2(,2())(,([11++-++++'-'=++-+++++-=??++所以上式为+--+''=?++dx hx x x y e n nx x n n n )2()(11θdx x y x f h x y h x f h x y h x f n n n n x x n n n n))](,(2)(,2())2(,2([1++-++?+ 3218)(LMh h x y Lh e n n ≤+''≤+?中点公式的整体截断误差：dx y x f hy h x f x y x f y x y y x y n n x x n n n n n n n n)],(2,2())(,([)()(111?+++-+-=-++dxy x f hy h x f x y x f h x y h x f x y x f hx y h x f x y x f y x y n n n n n n n n x x n n n n n n n n))],(2,2()))(,(2)(,2()))(,(2)(,2())(,([)(1++-+++++-+-=?+因⽽n n n L h Lh R εεε)21(1+++≤+，R L h Lh n n +++≤-122)21(εε≤≤])21()21(1[2)21(1222222022-+++++++--+++n nL h Lh L h Lh Lh Lh RL h Lh ε )1(00-+≤--x X L x X L e LhR eε 5．略 6．略 7．略8．（1）欧拉法：2.0≤h ；四阶Runge-Kutta ⽅法：278.0≤h （2）欧拉法：3 54≤h ；四阶Runge-Kutta ⽅法：3556.5≤h（3）欧拉法：1≤h ；四阶Runge-Kutta ⽅法：278.0≤h 9．略 10．略习题21．略 2．略 3．略4．差分格式写成矩阵形式为：n n M n M n n n M n M n n e u u u u r t r r r t r r r t r r r t u u u u +?--------= --+-+-++12211221121212121 αβαααβαααβαααβ矩阵的特征值为：)cos(221Mj r r t j πααβλ+-?-=，要使格式稳定，则特征值须满⾜t c j ?+≤1λ，即21≤r α5．利⽤泰勒展式可以得到古典隐式差分格式的截断误差为)(2h t O +?。

实验4常微分方程数值解实验目的：1.练习数值积分的计算；2.掌握用MATLAB软件求微分方程初值问题数值解的方法；3.通过实例学习用微分方程模型解决简化的实际问题；4.了解欧拉方法和龙格——库塔方法的基本思想和计算公式，及稳定性等概念。

实验内容：3.小型火箭初始质量为1400kg，其中包括1080kg燃料，火箭竖直向上发射是燃料燃烧率为18kg/s，由此产生32000N的推力，火箭引擎在燃料用尽时关闭。

设火箭上升是空气阻力正比于速度的平方，比例系数为0.4kg/m，求引擎关闭瞬间火箭的高度，速度，加速度，及火箭到达最高点是的高度，速度和加速度，并画出高度，速度，加速度随时间变化的图形。

解答如下：这是一个典型的牛顿第二定律问题，分析火箭受力情况；先规定向上受力为正数建立数学模型：A燃料未燃尽前，在任意时刻（t<60s）火箭受到向上的-F=32000N，向下的重力G=mg，g=9.8，向下的阻力f=kv^2, k=0.4, v表示此时火箭速度；此时火箭收到的合力为F1=（F-mg-f）；火箭的初始质量为1400kg，燃料燃烧率为-18kg/s;此刻火箭质量为m=1400-18*t根据牛顿第二定律知，加速度a=F1/m=(F-mg-f)/(m-r*t)=(32000-(0.4.*v.^2)-9.8.*(1400-18.*t))由此可利用龙格-库塔方法来实现，程序实现如下Function [dx]=rocket[t,x] %建立名为rocket的方程m=1400;k=0.4;r=-18;g=9.8; %给出题目提供的常数值dx=[x(2);(32000-(k*x(2)^2)-g*(m+r*t))/(m+r*t)];%以向量的形式建立方程[a]=(32000-(k*x(2)^2)-g*(m+r*t))/(m+r*t); %给出a的表达式End;ts=0:60; %根据题目给定燃烧率计算出燃料燃尽的时间，确定终点x0=[0,0]; %输入x的初始值[t,x]=ode15s(@rocket,ts,x0); %调用ode15s计算[t,x];h=x(:,1);v=x(:,2);plot(t,x(:,1)),grid; %绘出火箭高度与时间的关系曲线title('h-t');xlabel('t/s');ylabel('h/m'),pause;plot(t,x(:,2)),grid ; %绘出火箭速度与时间的曲线关系title('v-t');xlabel('t/s');ylabel('v/m/s'),pause;a=(32000-(0.4.*v.^2)-9.8.*(1400-18.*t))/(1400-18.*t); plot(t,a),grid; %绘出火箭加速度与时间的曲线关系title('a-t');xlabel('t/s'),ylabel('a/m^2/s'),pause火箭高度随时间变化的曲线火箭速度随时间变化的曲线火箭加速度随时间变化的曲线数据过多，故截取部分如下第一列为时间，第二列为火箭高度，第三列为火箭速度由此可以，在t=60s时，即火箭燃料燃尽瞬间，引擎关闭瞬间，火箭将到达12912m的高度，速度为267,29m，加速度a=0.9m/s^2B燃料燃尽之后，与A 类似，分析受力如下火箭受到向上的F=0向下的重力G=mg，g=9.8，向下的阻力f=kv^2, k=0.4, v表示此时火箭速度；此时火箭收到的合力为F2=（-mg-f）；火箭的初始质量为320kg，恒定根据牛顿第二定律,加速度a=F2/m=-g-0.4v^2/320;程序实现如下function [ dx ] = rocket2( t,x ) %建立以rocket2为名的函数dx=[x(2);-9.8-0.4.*x(2).^2/320]; %以向量的形式建立方程ts=60:120; %给出初始时刻，估计终点时刻x0=[12190,267.26]; %给出x初始值[t,x]=ode15s(@rocket2,ts,x0); %调用ode15s计算[t,x]plot(t,x(:,1)),grid; %绘出火箭高度随时间变化的曲线title('h-t');xlabel('t/s'),ylabel('h/m'),pause;plot(t,x(:,2)),grid; %绘出火箭速度随时间的变化曲线title('v-t');xlabel('t/s'),ylabel('v/m/s'),pause;v=x(:,2);a=-9.8-0.4*v.^2/320; %给出加速度的具体表达式plot(t,a),grid; %绘出火箭加速度随时间变化的曲线title('a-t');xlabel('t/s'),ylabel('a/m^2/s'),pause得到的曲线图形如下火箭高度随时间的变化曲线从图中可以大致看出，最高点在13km左右，火箭速度随时间的变化曲线加速度随时间变化曲线如下数据表格大致如下从图表中可以看出，在71s左右速度到达0，即此时到达最高处，高度为13117m加速度为-9.8m/m/s^2;本题总结：这道题是典型的物理牛顿力学的题目，通过受力的正确分析，可以知道，以[h,v]为向量建立微分方程即可求解，h的微分是速度v，速度v的微分是加速度a解题过程中存在的难点是：取值步长不太容易确定，而且是哪种算法不确定，先用ode15s 速度较快，ode23s速度差不太多，其他两种速度较慢，等待时间较长5．一只小船渡过宽为d 的河流，目标是起点A 正对着的另一岸B 点。

数值分析计算实习题答案数值分析计算实习题答案数值分析是一门研究如何利用计算机对数学问题进行近似求解的学科。

在数值分析的学习过程中，实习题是一种重要的学习方式，通过实践来巩固理论知识，并培养解决实际问题的能力。

本文将为大家提供一些数值分析计算实习题的答案，希望能够帮助大家更好地理解和掌握数值分析的相关知识。

一、插值与拟合1. 已知一组数据点，要求通过这些数据点构造一个一次插值多项式，并求出在某一特定点的函数值。

答案：首先，我们可以根据给定的数据点构造一个一次插值多项式。

假设给定的数据点为(x0, y0), (x1, y1)，我们可以构造一个一次多项式p(x) = a0 + a1x，其中a0和a1为待定系数。

根据插值条件，我们有p(x0) = y0，p(x1) = y1。

将这两个条件代入多项式中，可以得到一个方程组，通过求解这个方程组，我们就可以确定a0和a1的值。

最后，将求得的多项式代入到某一特定点，就可以得到该点的函数值。

2. 已知一组数据点，要求通过这些数据点进行最小二乘拟合，并求出拟合曲线的表达式。

答案：最小二乘拟合是一种通过最小化误差平方和来找到最佳拟合曲线的方法。

假设给定的数据点为(x0, y0), (x1, y1)，我们可以构造一个拟合曲线的表达式y =a0 + a1x + a2x^2 + ... + anx^n，其中a0, a1, ..., an为待定系数。

根据最小二乘拟合原理，我们需要最小化误差平方和E = Σ(yi - f(xi))^2，其中yi为实际数据点的y值，f(xi)为拟合曲线在xi处的函数值。

通过求解这个最小化问题，我们就可以确定拟合曲线的表达式。

二、数值积分1. 已知一个函数的表达式，要求通过数值积分的方法计算函数在某一区间上的定积分值。

答案：数值积分是一种通过将定积分转化为数值求和来近似计算的方法。

假设给定的函数表达式为f(x)，我们可以将定积分∫[a, b]f(x)dx近似为Σwi * f(xi)，其中wi为权重系数，xi为待定节点。

微分方程数值解法实验报告姓名： 班级： 学号：一：问题描述求解边值问题：()2(sin cos cos sin (0,1)(0,1)0,(,)x y u e x y x y G u x y G ππππππ+⎧⎫∆=+⎪⎪∈=⨯⎨⎬⎪⎪=∈∂⎩⎭（x,y) 其精确解为)sin()sin(),()(y x e y x u y x πππ+=问题一：取步长h=k=1/64,1/128,作五点差分格式，用Jacobi 迭代法，Gauss_Seidel 迭代法，SOR 迭代法（w=1.45）。

求解差分方程，以前后两次重合到小数点后四位的迭代值作为解的近似值，比较三种解法的迭代次数以及差分解)128/1,64/1)(,(=h y x u h 与精确解的精度。

问题二：取步长h=k=1/64,1/128,作五点差分格式，用单参数和双参数PR 法解差分方程，近似到小数点后四位。

与SOR 法比较精度和迭代步数。

问题三：取步长h=k=1/64,1/128,作五点差分格式，用共轭梯度法和预处理共轭梯度法解差分方程，近似到小数点后四位。

与SOR法与PR 法比较精度和迭代步数。

二．实验目的：分别使用五点差分法（Jacobi 迭代，Gauss_Seidel 迭代，SOR迭代），PR 交替隐式差分法（单参数，双参数），共轭梯度法，预共轭梯度法分别求椭圆方程的数值解。

三．实验原理：(1) Jacobi 迭代法设线性方程组(1)的系数矩阵A 可逆且主对角元素均不为零,令 并将A 分解成(2) 从而(1)可写成令其中. (3) 以为迭代矩阵的迭代法(公式)(4)称为雅可比(Jacobi)迭代法(公式),用向量的分量来表示,(4)为(5) 其中为初始向量. (2) Guass-Seidel 迭代法由雅可比迭代公式可知,在迭代的每一步计算过程中是用的全部分量来计算的所有分量,显然在计算第i 个分量时,已经b Ax =nn a ,...,a ,a 2211()nn a ,...,a ,a diag D 2211=()D D A A +-=()b x A D Dx +-=11f x B x +=b D f ,A D I B 1111--=-=1B ()()111f x B x k k +=+⎩⎨⎧[],...,,k ,n ,...,i x a b a x n ij j )k (j j i i ii )k (i 21021111==∑-=≠=+()()()()()Tn x ,...x ,x x 002010=()k x ()1+k x ()1+k i x计算出的最新分量没有被利用，从直观上看,最新计算出的分量可能比旧的分量要好些.因此，对这些最新计算出来的第次近似的分量加以利用,就得到所谓解方程组的高斯—塞德（Gauss-Seidel ）迭代法.把矩阵A 分解成(6)其中,分别为的主对角元除外的下三角和上三角部分,于是,方程组(1)便可以写成 即其中(7)以为迭代矩阵构成的迭代法(公式)(8)称为高斯—塞德尔迭代法(公式),用 量表示的形式为(3) SOR 迭代(4) 交替方向迭代法（PR 法）迭代格式为：()()1111+-+k i k x ,...,x 1+k ()1+k x ()1+k j x U L D A --=()nn a ,...,a ,a diag D 2211=U ,L --A ()b Ux x L D +=-22f x B x +=()()b L D f ,U L D B 1212---=-=2B ()()221f x B x k k +=+⎩⎨⎧[],...,,k ,n ,,i x a x a b a x i j n i j )k (j ij )k (j ij i ii )k (i 21021111111==∑∑--=-=+=++Λ))1(()(1D R L D T ωωω-+-=-b )(1--=L D d ωωhu πμωcos )11/(22opt =-+=2121,,1,1,1,,122L L L L u u u L u u u j i j i j i j i j i j i +==+-=+-+-+-对于单参数PR 法，对于多参数，(5) 共轭梯度法 算法步骤如下： ［预置步］任意，计算，并令取：指定算法终止常数，置，进入主步；［主步］ （１）如果，终止算法，输出；否则下行；（２）计算：（３）计算：（４）置，转入（１）．(6) 预共轭梯度法b uL I uL I b u L I uL I k k k k k k k k k k ττττττ+-=++-=++++211122211)()()()(hh optπτsin 22=2sin a ....2,1)11(421k 221h k a h k πρρτ==+-=--其中［预置步］任意，计算，并令取：指定算法终止常数，置，进入主步；［主步］（１）计算：，（２）如果，转入（3）．否则，终止算法，输出计算结果（３）计算：（４）置，转入（1）注：在算法［主步］中，引入变量，及，可以简化计算。

计算实验课微分方程数值解法数值计算实验题目一、常微分方程部分：1．使用四阶Runge-Kutta 方法求解如下初值问题的近似解，并将结果与实际值进行比较。

2．使用四阶Adams 预估校正算法（PECP 和PMECME 方案），初始值用四阶Runge-Kutta 方法提供，并将结果与实际值进行比较。

()21u t u -+='，32≤≤t ，()12=u ；精度510-=ε，5.0=h 。

实际解11u t t=+-。

tuu +='1，21≤≤t ，()21=u ；精度510-=ε，2.0=h 。

实际解2ln +=t t u 。

二、偏微分方程部分：1．用有限差分法求解如下Poisson 方程(),cos3sin ,u x y x y π-∆=，0x π<<，10<<y边界条件为： ()(),0,10,u x u x ==01x ≤≤； ()()0,,0,x x u y u y π==10≤≤y 取1,h Nπ=和21,h N=作矩形剖分，网格节点为1i x ih =，2j y jh =，i ,j =0,1,…,N 。

差分格式为1,,1,,1,,11222cos3sin i j i j i j i j i j i j i j u u u u u u x y h h π+-+--+-+⎛⎫-+= ⎪⎝⎭,1,2,,1i j N =-L 边界条件为： 00,0,,,i iN u u i N ===L01,1,,1,j j u u j N ==-L 1,1,,1,N j N j u u j N -==-L 结果与精确解()()12,9cos3sin u x y x y ππ-=+进行比较。

求解方案：依次令4,8,16,32N =，取6位小数计算。

用消元法求解，并就(),,44j i j x y π⎛⎫= ⎪⎝⎭，,1,2,3i j =处列出差分解与精确解。

其次，就N =32，0.25，0.5，0.75及i =0，2，4，…，30，32画出差分解曲线。