微分方程数值解实验3
微分方程的解
微分方程数值解
dy f ( x, y ) x [ a , b ] dx y ( a ) y0
ode23tb 刚性 刚性 分；精度中等 尝试使用 单步法；2 阶Rosebrock 算 当精度较低时，计算时 间比 ode15s 短 法；低精度 梯形算法；低精度 当精度较低时，计算时 间比ode15s短
注意: 1、在解n个未知函数的方程组时，y0和y均为n维向量， m-文件中的待解方程组应以x的分量形式写成. 2、如果需求解的问题是高阶常微分方程，则需将其化 为一阶常微分方程组，此时需用函数文件来定义该常微 分方程组. 3、odefun 为显式常微分方程，可以用命令inline 定义，或在函数文件中定义，然后通过函数句柄调用.
k k
几何意义
亦称为欧拉折线法
/* Euler’s polygonal arc method*/
Euler 折线法举例
例6：用 Euler 法解初值问题 2x dy y 2 y x [0, 2 ] dx y(0) 1 解：取步长 h = (2 - 0)/n = 2/n，得差分方程
n
等距剖分： a x 0 x1 x 2 x n 1 x n b 步长：h x k 1 x k ( b a ) / n , 差商代替微商
y ( x k 1 ) y ( x k ) dy y ( x k 1 ) y ( x k ) h y '( x k ) dx x h k y0 y ( x 0 ) k = 0, 1, 2, ..., n-1 得方程组： y k 1 y k h f ( x k , y k ) x x h y 是 y (x ) 的近似 k k 1
数值计算实验报告-欧拉法常微分方程
数学与计算科学学院实验报告实验项目名称欧拉法解常微分方程所属课程名称数值计算实验类型验证型实验日期2012-6- 4班级隧道1002班学号201008020233姓名李彬彬成绩一、实验概述：【实验目的】 通过运用相关的数值计算软件，解决最基本的常微分方程的数值计算，并且能够熟练的运用这种方法。
【实验原理】 欧拉法1．对常微分方程初始问题(9.2))((9.1)),(00⎪⎩⎪⎨⎧==y x y y x f dxdy用数值方法求解时，我们总是认为(9.1)、(9.2)的解存在且唯一。
欧拉法是解初值问题的最简单的数值方法。
从（9.2）式由于y (x 0) = y 0已给定，因而可以算出),()('000y x f x y =设x 1 = h 充分小，则近似地有：),()(')()(00001y x f x y hx y x y =≈-(9.3)记 ,n ,,i x y y i i 10 )(== 从而我们可以取),(0001y x hf y y ==作为y (x 1)的近似值。
利用y 1及f (x 1, y 1)又可以算出y (x 2)的近似值：),(1112y x hf y y +=一般地，在任意点x n +1 = (n + 1)h 处y (x )的近似值由下式给出),(1n n n n y x hf y y +=+（9.4）这就是欧拉法的计算公式，h 称为步长。
不难看出，近似解的误差首先是由差商近似代替微商（见(9.3)）引起的，这种近似代替所产生的误差称为截断误差。
还有一种误差称为舍入误差，这种误差是由于利用（9.4）进行计算时数值舍入引起的。
【实验环境】Windows XP 环境下运行 NumericalAnalyse 软件二、实验内容：【实验方案】在区间[0,1]上以h=0.1为步长，分别用欧拉法与预估-校正法求初值问题y’=y-2x/y且 y|x=0 =1的数值解。
将上述方程输入到软件NumericalAnalyse中步骤如图选择常微分方程的数值解法。
《微分方程数值解》实验指导书2010
第1页共18页微分方程数值解实验指导书李光云数学与计算科学学院O 一0十月第2页共18页一、概述本课程实验指导书是根据陆金甫，关治编著的《偏微分方程数值解法（第二版）》编写的。
通过上机实验，可帮助学生理解偏微分方程数值计算方法的基本思想，巩固学生所学过的算法，同时让学生进一步学习和掌握Matlab软件在数值计算方面中的应用。
二、实验环境本书选择的实验环境是计算机以及软件Matlab（版本6.5以上）一套。
三、实验课时安排共7次实验，16课时，试验五为综合设计性试验4课时，其它实验2个课时。
四、实验要求上机完成实验指导书中所规定的内容，自行按实验指导书要求完成程序设计和调试，并提交每次实验的实验报告，附带算法程序清单和算法输出结果。
五、实验考核要求上机完成试验内容，并提交一份算法程序清单和数值结果。
第3页共18页实验一 双曲型方程的迎风格式、实验目的1、 掌握双曲型方程的迎风格式的算法思想 ，格式稳定的条件;2、 培养编程与上机调试能力。
、实验课时：2个课时 三、实验准备1、 了解偏微分方程的分类，以及双曲型方程的特点；2、 熟悉求解偏微分方程有限差分法的基本原理，双曲型方程迎风格式的推导和 特点，熟悉用迎风格式求解双曲型方程的算法步骤；3、 熟悉Matlab 软件的基本命令及其数值计算中的基本命令和函数。
4、 了解双曲型方程迎风格式稳定的条件。
四、实验内容用双曲型方程的迎风格式n nU j -U j 」0, hn 1 nnnU j —U jU j1—U j c a 0,Th其中，• ,h 分别为时间步长和空间步长。
求解初值问题賦”,ts,I u(x,0 )=u°(x ),其中,p!,x",U。
八 0,x 0.五、基本思想及主要步骤考虑对流方程:U :ua 0,x . ,t 0* *-t:x其a 为给定常数。
其迎风格式的基本思想是在方程中关于空间的偏导数用在特征线n 1 nU j Uja 0 a<0.计算〔0,1〕区间的u 值至t n 二0.5,并画出解的图形。
数学实验——常微分方程数值解
实验4 常微分方程数值解分1 黄浩 43一、实验目的1.掌握用MATLAB软件求微分方程初值问题数值解的方法；2.通过实例学习用微分方程模型解决简化的实际问题；3.了解欧拉方法和龙格-库塔方法的基本思想和计算公式，及稳定性等概念。
二、实验内容1.《数学实验》第一版（问题2）问题叙述：小型火箭初始重量为1400kg，其中包括1080kg燃料。
火箭竖直向上发射时燃料燃烧率为18kg/s，由此产生32000N的推力，火箭引擎在燃烧用尽时关闭。
设火箭上升时空气阻力正比于速度的平方，比例系数为0.4kg/m，求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。
模型转换及实验过程：（一）从发射到引擎关闭设火箭总质量为m，上升高度为h，瞬时速度为v，瞬时加速度为a，由燃料燃烧时间t=60s，可列如下的方程组：t∈[0,60]−9.8因此，上述方程为二元常微分方程组，选择t为自变量，h和v为因变量进行分析。
初值条件：h(0)=0 ,v(0)=0对上述模型，使用ode45()函数求数值解（程序见四.1、四.2），结果如下：由上表可知，引擎关闭瞬间，火箭的高度为12189.78m，速度为267.26m/s，加速度为0.9170m/s2，火箭至此已飞行60s而高度、速度、加速度随时间的变化曲线如下：(二)从引擎关闭到最高点设引擎关闭时，t1=0，由上一问的结果可知，h(0)=12189.78m，v(0)= 267.26m/s，m=320kg，则可列二元常微分方程组如下：9.8因此，可选择t1为自变量，h、v为因变量进行分析（程序见四.3、四.4），实验结果如下：由上表可知，当t1∈[11,12]时，v(t1)有零点，即该区间内某时刻火箭达到最高点。
再进行更细致的实验（程序略），设步长为0.01，观察该区间内v(t1)的零点，如下表所示：可以看出，当t1=11.30s，即总时间t=71.30s时，火箭达到最高点，高度为13115.36m，加速度为-9.8m/s2。
微分方程的数值解法
微分方程的数值解法微分方程是描述自然界中众多现象和规律的重要数学工具。
然而，许多微分方程是很难或者无法直接求解的，因此需要使用数值解法来近似求解。
本文将介绍几种常见的微分方程数值解法。
1. 欧拉方法欧拉方法是最简单的数值解法之一。
它将微分方程转化为差分方程，通过计算离散点上的导数来逼近原方程的解。
欧拉方法的基本思想是利用当前点的导数值来估计下一个点的函数值。
具体步骤如下：首先，将自变量区间等分为一系列的小区间。
然后，根据微分方程的初始条件，在起始点确定初始函数值。
接下来，根据导数的定义，计算每个小区间上函数值的斜率。
最后，根据初始函数值和斜率，递推计算得到每个小区间上的函数值。
2. 龙格-库塔方法龙格-库塔方法是一种常用的高阶精度数值解法。
它通过进行多次逼近和修正来提高近似解的准确性。
相比于欧拉方法，龙格-库塔方法在同样的步长下可以获得更精确的解。
具体步骤如下：首先，确定在每个小区间上的步长。
然后，根据微分方程的初始条件，在起始点确定初始函数值。
接下来，根据当前点的导数值，使用权重系数计算多个中间点的函数值。
最后，根据所有中间点的函数值，计算出当前点的函数值。
3. 改进欧拉方法（改进的欧拉-克罗默法）改进欧拉方法是一种中阶精度数值解法，介于欧拉方法和龙格-库塔方法之间。
它通过使用两公式递推来提高精度，并减少计算量。
改进欧拉方法相对于欧拉方法而言，增加了一个估计项，从而减小了局部截断误差。
具体步骤如下：首先，确定在每个小区间上的步长。
然后，根据微分方程的初始条件，在起始点确定初始函数值。
接下来，利用欧拉方法计算出中间点的函数值。
最后，利用中间点的函数值和斜率，计算出当前点的函数值。
总结：微分方程的数值解法为我们研究和解决实际问题提供了有力的工具。
本文介绍了欧拉方法、龙格-库塔方法和改进欧拉方法这几种常见的数值解法。
选择合适的数值解法取决于微分方程的性质以及对解的精确性要求。
在实际应用中，我们应该根据具体情况选择最合适的数值解法，并注意控制步长以尽可能减小误差。
微分方程数值解法实验报告
微分方程数值解法实验报告2篇微分方程数值解法实验报告（一）在实际科学与工程问题中，我们经常会遇到微分方程的求解。
然而，许多微分方程往往没有解析解，这就需要我们利用数值方法来获得近似解。
本篇实验报告将介绍两种常见的微分方程数值解法：欧拉方法和改进的欧拉方法。
一、欧拉方法欧拉方法是最简单的微分方程数值解法之一。
其基本原理为离散化微分方程，将微分方程中的导数用差商代替。
设要求解的微分方程为dy/dx = f(x, y)，步长为h，则可用以下公式进行递推计算：y_{n+1} = y_n +hf(x_n, y_n)二、算法实现为了对欧拉方法进行数值实验，我们以一阶线性常微分方程为例：dy/dx = x - y, y(0) = 1步骤如下：（1）选择合适的步长h和求解区间[a, b]，这里我们取h=0.1，[a, b] = [0, 1]；（2）初始化y_0 = 1；（3）利用欧拉方法递推计算y_{n+1} = y_n + 0.1(x_n - y_n)；（4）重复步骤（3），直到x_n = 1。
三、实验结果与讨论我们通过上述步骤得到了在[0, 1]上的近似解y(x)。
下图展示了欧拉方法求解的结果。
从图中可以看出，欧拉方法得到的近似解与精确解有一定的偏差。
这是因为欧拉方法只是通过递推计算得到的近似解，并没有考虑到更高阶的误差。
所以在需要高精度解时，欧拉方法并不理想。
四、改进的欧拉方法针对欧拉方法的不足，我们可以考虑使用改进的欧拉方法（也称为改进的欧拉-柯西方法）。
它是通过利用前后两个步长欧拉方法得到的结果作为差商的中间项，从而提高了求解精度。
一阶线性常微分方程的改进欧拉方法可以表示为：y_{n+1} = y_n + h(f(x_n, y_n) + f(x_{n+1}, y_n + hf(x_n,y_n)))/2五、算法实现与结果展示将改进的欧拉方法应用于前述的一阶线性常微分方程，我们同样选择h=0.1，[a, b] = [0, 1]。
微分方程数值解实验报告
微分方程数值解实验报告实验目的：掌握微分方程数值解的基本方法，能够利用计算机编程求解微分方程。
实验原理：微分方程是自然科学与工程技术中常见的数学模型，它描述了变量之间的关系及其随时间、空间的变化规律。
解微分方程是研究和应用微分方程的基础，但有很多微分方程无法找到解析解，只能通过数值方法进行求解。
本实验采用欧拉方法和改进的欧拉方法求解微分方程的初值问题：$$\begin{cases}\frac{dy}{dt}=f(t,y)\\y(t_0)=y_0\end{cases}$$其中，$f(t,y)$是给定的函数，$y(t_0)=y_0$是已知的初值条件。
欧拉方法是最基本的数值解法，其步骤如下：1.给定$t_0$和$y_0$2.计算$t_{i+1}=t_i+h$，其中$h$是步长3. 计算$y_{i+1}=y_i+hf(t_i,y_i)$4.重复步骤2、3直到达到终止条件改进的欧拉方法是对欧拉方法进行改进，通过利用函数$y(t)$在$t+\frac{1}{2}h$处的斜率来更准确地估计$y_{i+1}$，其步骤如下：1.给定$t_0$和$y_0$2.计算$t_{i+1}=t_i+h$，其中$h$是步长3. 计算$y_*=y_i+\frac{1}{2}hf(t_i,y_i)$4. 计算$y_{i+1}=y_i+hf(t_i+\frac{1}{2}h,y_*)$5.重复步骤2、3、4直到达到终止条件实验步骤：1.编写程序实现欧拉方法和改进的欧拉方法2.给定微分方程和初值条件3.设置步长和终止条件4.利用欧拉方法和改进的欧拉方法求解微分方程5.比较不同步长下的数值解与解析解的误差6.绘制误差-步长曲线，分析数值解的精度和收敛性实验结果：以一阶常微分方程$y'=3ty+t$为例，给定初值$y(0)=1$，取步长$h=0.1$进行数值求解。
利用欧拉方法求解微分方程得到的数值解如下：\begin{array}{cccc}t & y_{\text{exact}} & y_{\text{Euler}} & \text{误差} \\ \hline0.0&1.000&1.000&0.000\\0.1&1.035&1.030&0.005\\0.2&1.104&1.108&0.004\\0.3&1.212&1.217&0.005\\0.4&1.360&1.364&0.004\\0.5&1.554&1.559&0.005\\0.6&1.805&1.810&0.005\\0.7&2.131&2.136&0.005\\0.8&2.554&2.560&0.006\\0.9&3.102&3.107&0.006\\1.0&3.807&3.812&0.005\\\end{array}利用改进的欧拉方法求解微分方程得到的数值解如下：\begin{array}{cccc}t & y_{\text{exact}} & y_{\text{Improved Euler}} & \text{误差} \\\hline0.0&1.000&1.000&0.000\\0.1&1.035&1.035&0.000\\0.2&1.104&1.103&0.001\\0.3&1.212&1.211&0.001\\0.4&1.360&1.358&0.002\\0.5&1.554&1.552&0.002\\0.6&1.805&1.802&0.003\\0.7&2.131&2.126&0.005\\0.8&2.554&2.545&0.009\\0.9&3.102&3.086&0.015\\1.0&3.807&3.774&0.032\\\end{array}误差-步长曲线如下：实验结论：通过对比欧拉方法和改进的欧拉方法的数值解与解析解的误差，可以发现改进的欧拉方法具有更高的精度和收敛性。
河北工业大学数值分析实验三实验四实验报告
数值分析实验报告指导老师：宛艳萍姓名：班级：学号：实验三 复化辛卜生法，龙贝格法1.实验名称：复化辛卜生法，龙贝格法2.实验目的1）通过实际计算体会各种方法的精确度。
2）会编写用复化辛卜生、龙贝格算法求定积分的程序。
3.算法描述1）用复化辛卜生法计算积分 dxx I ⎰+=12)1/(1算法：复化辛卜生公式为S n =h/6∑∑+-=+++)]()2/(4)([11k k kn k x f h x f xf ,计算过程为：1．令,/)(n a b h -= ),2/(1h a f s +=;02=s2．对1,,2,1-=n k计算),2/(11h kh a f s s +++=)(22kh a f s s ++=3．))(24)((6/21b f s s a f h s +++= 。
2）龙贝格算法计算dxxI ⎰+=102)1/(156e ε=-算法)((12/12∑-=++=n k k n n n x f h T T ;/)(n a b h n -= n k h k x )2/1(2/1+=+)(3/122n n n n T T T S -+= )_(15/122n n n n S S S C +=)(63/122n n n n C C C R -+=用事后估计法控制精度2|5e -6n n R R -< 。
4.源程序：1）/* 用复化辛卜生公式求积分 */ #include "stdio.h" float fx(float x){double f;f=1.0/(1.0+x*x); return f; } double fs(int n){double a=0.0,b=1.0,h,s,s1,s2=0; int i;h=(b-a)/n; s1=fx(a+h/2); for(i=1;i<n;i++){s1=s1+fx(a+i*h+h/2); s2=s2+fx(a+i*h);}s=(h/6.0)*(fx(a)+fx(b)+4*s1+2*s2);return s;}void main(){printf("实验三复化辛卜生法计算机112 耿向飞学号：112434\n");printf("s(2)=%lf\ns(4)=%lf\ns(8)= %lf",fs(2),fs(4),fs(8));}2）/* 龙贝格法 */#include "stdio.h"#include "math.h"#define E 2.71828182//被积函数f(x)double fx(double x){double f;f=1/(1+x*x);return f;}//梯形公式求tndouble tx(int n){double s3=0.0,h,t,b=1.0,a=0.0;int i;h=(b-a)/n;for(i=1;i<n;i++)s3=s3+fx(i*h);t=(h/2)*(fx(a)+fx(b)+2*s3);return t;} double s(int n){double s;s=tx(2*n)+(1.0/3.0)*(tx(2*n)-tx(n ));return s;}double c(int n){double c;c=s(2*n)+(1.0/15.0)*(s(2*n)-s(n)) ;return c;}double r(int n){double r;r=c(2*n)+(1.0/63.0)*(c(2*n)-c(n)) ;return r;}void main(){double rr,pp;int n=1;rr=r(n);pp=r(2*n)-r(n);printf("实验三龙贝格法计算机112 耿向飞学号：112434\n");printf("结果为：%.15lf 误差小于等于： %.15lf",rr,pp);}5.运行结果1）复化辛卜生公式2）龙贝格算法6.对算法的理解与分析：复化辛卜生公式和龙贝格算法适用于求数值积分，而且都能提高计算积分的精度龙贝格算法其实是在复化辛卜生公式递推的基础之上生成的一种精度高，而且收敛速度也较快的一种算法。
实验报告七常微分方程初值问题的数值解法
浙江大学城市学院实验报告课程名称数值计算方法实验项目名称常微分方程初值问题的数值解法 实验成绩指导老师签名日期2015/12/16 一.实验目的和要求1． 用Matlab 软件掌握求微分方程数值解的欧拉方法和龙格－库塔方法； 2． 通过实例学习用微分方程模型解决简化的实际问题;二.实验内容和原理编程题2-1要求写出Matlab 源程序m 文件,并有适当的注释语句；分析应用题2-2,2-3,2-4,2-5要求将问题的分析过程、Matlab 源程序和运行结果和结果的解释、算法的分析写在实验报告上; 2-1 编程编写用向前欧拉公式和改进欧拉公式求微分方程数值解的Matlab 程序,问题如下：在区间[],a b 内(1)N +个等距点处,逼近下列初值问题的解,并对程序的每一句添上注释语句; Euler 法y=eulera,b,n,y0,f,f1,b1改进Euler 法y=eulerproa,b,n,y0,f,f1,b1 2-2 分析应用题假设等分区间数100n =,用欧拉法和改进欧拉法在区间[0,10]t ∈内求解初值问题()()20(0)10y t y t y '=-⎧⎨=⎩并作出解的曲线图形,同时将方程的解析解也画在同一张图上,并作比较,分析这两种方法的精度; 2-3 分析应用题用以下三种不同的方法求下述微分方程的数值解,取10h = 画出解的图形,与精确值比较并进行分析; 1欧拉法； 2改进欧拉法； 3龙格－库塔方法；2-4 分析应用题考虑一个涉及到社会上与众不同的人的繁衍问题模型;假设在时刻t 单位为年,社会上有人口()x t 人,又假设所有与众不同的人与别的与众不同的人结婚后所生后代也是与众不同的人;而固定比例为r 的所有其他的后代也是与众不同的人;如果对所有人来说出生率假定为常数b ,又如果普通的人和与众不同的人的婚配是任意的,则此问题可以用微分方程表示为：其中变量()()()i p t x t x t =表示在时刻t 社会上与众不同的人的比例,()i x t 表示在时刻t 人口中与众不同的人的数量;1假定(0)0.01,0.02p b ==和0.1r =,当步长为1h =年时,求从0t =到50t =解()p t 的近似值,并作出近似解的曲线图形;2精确求出微分方程的解()p t ,并将你当50t =时在分题b 中得到的结果与此时的精确值进行比较; MATLAB 相关函数求微分方程的解析解及其数值的代入dsolve‘egn1’,‘egn2’,‘x ’ subsexpr,{x,y,…},{x1,y1,…}其中‘egn i ’表示第i 个方程,‘x ’表示微分方程中的自变量,默认时自变量为t ; subs 命令中的expr 、x 、y 为符合型表达式,x 、y 分别用数值x1、x2代入; >>symsxyz>>subs'x+y+z',{x,y,z},{1,2,3} ans= 6>>symsx>>subs'x^2',x,2 ans= 4>>s=dsolve‘12Dy y ∧=+’,‘(0)1y =’,‘x ’ ans= >>symsx >>subss,x,2 ans=右端函数(,)f x y 的自动生成f=inline ‘expr ’,’var1’,‘var2’,……其中’expr ’表示函数的表达式,’var1’,‘var2’表示函数表达式中的变量,运行该函数,生成一个新的函数表达式为fvar1,var2,……; >>f=inline'x+3y','x','y' f=Inlinefunction: fx,y=x+3y >>f2,3 ans= 114,5阶龙格－库塔方法求解微分方程数值解t,x=ode45f,ts,x0,options其中f 是由待解方程写成的m 文件名；x0为函数的初值；t,x 分别为输出的自变量和函数值列向量,t的步长是程序根据误差限自动选定的;若ts=t0,t1,t2,…,tf,则输出在自变量指定值,等步长时用ts=t0:k:tf,输出在等分点；options 用于设定误差限可以缺省,缺省时设定为相对误差310-,绝对误差610-,程序为：options=odeset ‘reltol ’,rt,’abstol ’,at,这里rt,at 分别为设定的相对误差和绝对误差;常用选项见下表;选项名 功能 可选值 省缺值 AbsTol 设定绝对误差正数 RelTol 设定相对误差 正数InitialStep 设定初始步长 正数 自动 MaxStep设定步长上界正数MaxOrder 设定ode15s 的最高阶数 1,2,3,4,5 5 Stats 显示计算成本统计 on,off off BDF 设定ode15s 是否用反向差分on,offoff例：在命令窗口执行>>odefun =inline ‘2*y t y -’,‘t ’,‘y ’；>>[],45(,[0,4],1)t y ode odefun =；ans=>>t y ‘o-’,%解函数图形表示>>45(,[0,4],1)ode odefun %不用输出变量,则直接输出图形 >>[],45(,0:4,1)t y ode odefun =；[],t yans=三.操作方法与实验步骤包括实验数据记录和处理2-1编程编写用向前欧拉公式和改进欧拉公式求微分方程数值解的Matlab 程序,问题如下：在区间[],a b 内(1)N +个等距点处,逼近下列初值问题的解,并对程序的每一句添上注释语句; Euler 法y=eulera,b,n,y0,f,f1,b1改进Euler 法y=eulerproa,b,n,y0,f,f1,b1Euler 法y=eulera,b,n,y0,f,f1,b1 y=zeros1,n+1; y1=y0; h=b-a/n; x=a:h:b; fori=1:n; yi+1=yi+hfxi,yi; end plotx,y holdon%求微分方程的精确解 x1=linspacea,b,100; '精确解为' s=dsolvef1,b1,'x' symsxy1=zeros1,100; for i=1:100y1i=subss,x,x1i; endplotx1,y1,'r'title'红色代表精确解'改进Euler 法y=eulerproa,b,n,y0,f,f1,b1 %求微分方程的数值解 y=zeros1,n+1; y1=y0; h=b-a/n; x=a:h:b; fori=1:n; T1=fxi,yi; T2=fxi+1,yi+hT1; yi+1=yi+h/2T1+T2; end plotx,y holdon%求微分方程的精确解 x1=linspacea,b,100; '精确解为' s=dsolvef1,b1,'x' symsxy1=zeros1,100; fori=1:100 y1i=subss,x,x1i; endplotx1,y1,'r'title'红色代表精确解' 2-2分析应用题假设等分区间数100n =,用欧拉法和改进欧拉法在区间[0,10]t ∈内求解初值问题()()20(0)10y t y t y '=-⎧⎨=⎩并作出解的曲线图形,同时将方程的解析解也画在同一张图上,并作比较,分析这两种方法的精度;1向前欧拉法>>euler0,10,100,10,inline'y-20','x','y','Dy=y-20','y0=10' ans= 精确解为 s= 20-10expx ans= +005Columns1through8（2）改进欧拉法>>eulerpro0,10,100,10,inline'y-20','x','y','Dy=y-20','y0=10' ans= 精确解为 s= 20-10expx ans= +005Columns1through8改进欧拉法的精度比向前欧拉法更高; 2-3分析应用题用以下三种不同的方法求下述微分方程的数值解,取10h = 画出解的图形,与精确值比较并进行分析; 1欧拉法； 2改进欧拉法；2-4分析应用题考虑一个涉及到社会上与众不同的人的繁衍问题模型;假设在时刻t 单位为年,社会上有人口()x t 人,又假设所有与众不同的人与别的与众不同的人结婚后所生后代也是与众不同的人;而固定比例为r 的所有其他的后代也是与众不同的人;如果对所有人来说出生率假定为常数b ,又如果普通的人和与众不同的人的婚配是任意的,则此问题可以用微分方程表示为：其中变量()()()i p t x t x t =表示在时刻t 社会上与众不同的人的比例,()i x t 表示在时刻t 人口中与众不同的人的数量;1假定(0)0.01,0.02p b ==和0.1r =,当步长为1h =年时,求从0t =到50t =解()p t 的近似值,并作出近似解的曲线图形;2精确求出微分方程的解()p t ,并将你当50t =时在分题b 中得到的结果与此时的精确值进行比较;1>>euler0,50,50,,inline'','t','p','Dp=','p0= 1' ans= 精确解为 s=1-99/100expx/500 ans=Columns1through82>>dsolve'Dp=','p0=','t' ans=1-99/100expt/500 >>1-99/100exp ans=与欧拉法求得的精确值差0,0001四.实验结果与分析。
偏微分方程数值解上机实验报告(matlab做的)
偏微分方程数值解法上机报告（一）一、实验题目：用Ritz-Galerkin 方法求解边值问题2u '',01(0)0,(1)1u x x u u ⎧-+=<<⎨==⎩的第n 次近似()n u x ，基函数()sin(),1,2,...,i x i x i n ϕπ==.二、实验目的：通过本次上机实验，理解求解初值问题的变分问题的最重要的近似解法——Ritz-Galerkin 方法，以便为学习有限元法打好基础。
此外，要熟悉用Matlab 解决数学问题的基本编程方法，提高运用计算机解决问题的能力。
三、实验代码：n=5;syms x;for i=1:np(i)=sin(i*pi*x);q(i)=-i^2*pi^2*sin(i*pi*x);endfor i=1:nb(i)=2*int(p(i),0,1);for j=1:nA(i,j)=int((-q(j)+p(j))*p(i),0,1);endendt=inv(A)*b'四、运行结果：t=2251799813685248/3059521645650671/pi281474976710656/9481460623939047/pi281474976710656/43582901062631895/pi五、总结：通过本次上机，我了解了Ritz-Galerkin 方程 n j j p f cj p i p a n i i ,...,2,1)),(,())(),((1==∑=，明白了用Ritz-Galerkin 方法解决边值问题的变分问题的基本原理，并接近一步提高自己的编程动手能力，受益匪浅。
偏微分方程数值解法上机报告（二）一、 实验题目：用线性元求下列边值问题的数值解2''2sin ,0142(0)0,'(1)0y y x x y y ππ⎧-+=<<⎪⎨⎪==⎩二、 实验目的：通过本次上机，熟悉和掌握用Galerkin 法观点出发导出的求解处置问题数值解的线性有限元法。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
-4.671818487622659e-012
步长
数值解
精确解
误差
0.2
3.71824930851150
3.718281828459046
-3.251994754860377e-005
0.1
3.71828096531189
3.718281828459046
-8.631471573039562e-007
1.790293186587633e-006
0.2
3.025527878319426
3.718281828459046
1.313731487551095e-005
function[y,errorPMCM]=PMECM(fun1,x0,y0,h,N)
x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;
《常微分方程数值解》实验报告
实验序号：3实验项目名称：线性多步方法之Adams方法
学 号
姓 名点
K7-405
指导教师
汪玉霞
时间
2014、4、24
一、实验目的及要求
1.让学生掌握用Adams内插方法
与外插方法 及Adams预估校正方
法求解常微分方程初值问题.
2.培养Matlab编程与上机调试能力.
x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;
forn=1:N
x(n+1)=x(n)+h;
ifn<4
k1=h*feval(fun1,x(n),y(n));
k2=h*feval(fun1,x(n)+0.5*h,y(n)+k1/2);
k3=h*feval(fun1,x(n)+0.5*h,y(n)+k2/2);
0.01
3.71828182845437
3.718281828459046
-4.671818487622659e-012
5、实验心得与实验结论
六、教师评语成绩
签名：
日期：
实验原理：用同阶精度的四步Adams外插方法求得的 作为预估值，再用三步Adams内插方法计算得到 的校正值，即利用了显式方法（外插方法）的便于计算的特点，又达到了隐式格式（内插方法）的计算精度
三、实验内容与步骤
实验内容：
1.取步长h=0.10用Adams三步内插方法和四步外插方法求解初值问题 并比较在同一点处的数值解的精度。（微分方程精确解为 ）或者用三步Adams内插方法和三步Adams外插方法分别计算同一点处得数值解，并比较精度。
ypm(n+1)=yp(n+1)+(251/270)*(y(n)-yp(n));
f(n+1)=h*feval(fun1,x(n+1),ypm(n+1));
yc(n+1)=y(n)+(9*f(n+1)+19*k1-5*k2+k3)/24;
y(n+1)=yc(n+1)-(19/270)*(yc(n+1)-yp(n+1));
y = 3.718282081879897
errorPC =2.534208514148872e-007
步长
数值解
精确解
误差（取绝对值）
0.05
3.718282081879897
3.718281828459046
2.534208514148872e-007
0.1
3.718283618752232
3.718281828459046
f(n+1)=h*feval(fun1,x(n+1),y(n+1));
end
end
errorPMCM=y(N+1)-1-exp(1);
>> [y,errorPMCM]=PMECM(fun1,0,1,0.2,5)
y =
3.71824930851150
errorPMCM =
-3.251994754860377e-005
yp=zeros(1,N+1);ypm=zeros(1,N+1);
yc=zeros(1,N+1);f=zeros(1,N+1);
forn=1:N
x(n+1)=x(n)+h;
ifn<4
k1=h*feval(fun1,x(n),y(n));
k2=h*feval(fun1,x(n)+0.5*h,y(n)+k1/2);
k3=h*feval(fun1,x(n)+0.5*h,y(n)+k2/2);
k4=h*feval(fun1,x(n)+h,y(n)+k3);
y(n+1)=y(n)+(k1+2*k2+2*k3+k4)/6;
yp(n+1)=y(n+1);ypm(n+1)=y(n+1);yc(n+1)=y(n+1);
else
errorPC =
-1.313731487551095e-005
> [y,errorPC]=Adams(fun1,0,1,0.1,10)
y =3.718283618752232
errorPC =1.790293186587633e-006
>> [y,errorPC]=Adams(fun1,0,1,0.05,20)
>> [y,errorPMCM]=PMECM(fun1,0,1,0.1,10)
y =
3.71828096531189
errorPMCM =
-8.631471573039562e-007
>>[y,errorPMCM]=PMECM(fun1,0,1,0.01,100)
y =
3.71828182845437
2 .用Adams预估校正修正格式求解
实验步骤
1．建立命令文件,保存为Adams.m
function [x,y]=Adams(fun1,x0,y0,h,N)
2．建立右端函数文件fun1.m
3.对具体的参数运行PMCM.m
4．对以上程序作局部修改，计算第1.2题
四、实验程序、结果与数据处理
function[y,errorPC]=Adams(fun1,x0,y0,h,N)
k4=h*feval(fun1,x(n)+h,y(n)+k3);
y(n+1)=y(n)+(k1+2*k2+2*k3+k4)/6;
else
k1=h*feval(fun1,x(n),y(n));
k2=h*feval(fun1,x(n-1),y(n-1));
k3=h*feval(fun1,x(n-2),y(n-2));
fun1=inline('-x+y+1','x','y');
>> [y,errorPC]=Adams(fun1,0,1,0.2,5)
y =1.000000000000000 1.421400000000000 1.891817960000000 2.422106456344000 3.025527878319426 3.718268691144170
k1=h*feval(fun1,x(n),y(n));f(n)=k1;
k2=h*feval(fun1,x(n-1),y(n-1));
k3=h*feval(fun1,x(n-2),y(n-2));
k4=h*feval(fun1,x(n-3),y(n-3));
yp(n+1)=y(n)+(55*k1-59*k2+37*k3-9*k4)/24;
3.理解并实现Adams预估校正修正格式：
二、实验设备要求与实验原理
1、支持Intel PentiumⅢ及其以上CPU，内存256MB以上、硬盘1GB以上容量的微机；软件配有Windows98/2000/XP操作系统及MATLAB软件；
2、实验过程中，务必将实验所得数据列表，分析数据，得出实验结论。（建议同时网上提交电子版实验报告：780668212@）
k4=h*feval(fun1,x(n-3),y(n-3));
yp=y(n)+(55*k1-59*k2+37*k3-9*k4)/24;
k0=h*feval(fun1,x(n+1),yp);
y(n+1)=y(n)+(9*k0+19*k1-5*k2+k3)/24;
end
end
errorPC=y(N+1)-1-exp(1);