常微分方程数值解实验报告
常微分方程数值解实验报告
常微分方程数值解实验报告学院：数学与信息科学专业：信息与计算科学姓名：郑思义学号：201216524课程：常微分方程数值解实验一：常微分方程的数值解法1、分别用Euler 法、改进的Euler 法（预报校正格式）和S —K 法求解初值问题。
（h=0.1）并与真解作比较。
⎩⎨⎧=++-=10(1y'）y x y 1.1实验代码：%欧拉法function [x,y]=naeuler(dyfun,xspan,y0,h)%dyfun 是常微分方程，xspan 是x 的取值范围，y0是初值，h 是步长 x=xspan(1):h:xspan(2); y(1)=y0;for n=1:length(x)-1y(n+1)=y(n)+h*feval(dyfun,x(n),y(n)); end%改进的欧拉法function [x,m,y]=naeuler2(dyfun,xspan,y0,h)%dyfun 是常微分方程，xspan 是x 的取值范围，y0是初值，h 是步长。
%返回值x 为x 取值，m 为预报解，y 为校正解 x=xspan(1):h:xspan(2); y(1)=y0;m=zeros(length(x)-1,1); for n=1:length(x)-1 k1=feval(dyfun,x(n),y(n)); y(n+1)=y(n)+h*k1; m(n)=y(n+1);k2=feval(dyfun,x(n+1),y(n+1)); y(n+1)=y(n)+h*(k1+k2)/2; end%四阶S —K 法function [x,y]=rk(dyfun,xspan,y0,h)%dyfun 是常微分方程，xspan 是x 的取值范围，y0是初值，h 是步长。
x=xspan(1):h:xspan(2); y(1)=y0;for n=1:length(x)-1 k1=feval(dyfun,x(n),y(n));k2=feval(dyfun,x(n)+h/2,y(n)+(h*k1)/2); k3=feval(dyfun,x(n)+h/2,y(n)+(h*k2)/2); k4=feval(dyfun,x(n)+h,y(n)+h*k3);y(n+1)=y(n)+(h/6)*(k1+2*k2+2*k3+k4);end%主程序x=[0:0.1:1];y=exp(-x)+x;dyfun=inline('-y+x+1');[x1,y1]=naeuler(dyfun,[0,1],1,0.1);[x2,m,y2]=naeuler2(dyfun,[0,1],1,0.1);[x3,y3]=rk(dyfun,[0,1],1,0.1);plot(x,y,'r',x1,y1,'+',x2,y2,'*',x3,y3,'o');xlabel('x');ylabel('y');legend('y为真解','y1为欧拉解','y2为改进欧拉解','y3为S—K解','Location','NorthWest');1.2实验结果：x 真解y 欧拉解y1 预报值m 校正值y2 S—K解y30.0 1.0000 1.0000 1.0000 1.00000.1 1.0048 1.0000 1.0000 1.0050 1.00480.2 1.0187 1.0100 1.0145 1.0190 1.01870.3 1.0408 1.0290 1.0371 1.0412 1.04080.4 1.0703 1.0561 1.0671 1.0708 1.07030.5 1.1065 1.0905 1.1037 1.1071 1.10650.6 1.1488 1.1314 1.1464 1.1494 1.14880.7 1.1966 1.1783 1.1945 1.1972 1.19660.8 1.2493 1.2305 1.2475 1.2500 1.24930.9 1.3066 1.2874 1.3050 1.3072 1.30661.0 1.3679 1.3487 1.3665 1.3685 1.36792、选取一种理论上收敛但是不稳定的算法对问题1进行计算，并与真解作比较。
实验报告——常微分方程的数值解法
实验报告实验项目名称常微分方程的数值解法实验室数学实验室所属课程名称微分方程数值解实验类型上机实验实验日期2013年3月11日班级10信息与计算科学学号2010119421姓名叶达伟成绩实验概述：【实验目的及要求】运用不同的数值解法来求解具体问题，并通过具体实例来分析比较各种常微分方程的数值解法的精度，为以后求解一般的常微分方程起到借鉴意义。
【实验原理】各种常微分方程的数值解法的原理，包括Euler法，改进Euler法，梯形法，Runge-Kutta方法，线性多步方法等。
【实验环境】（使用的软硬件）Matlab软件实验内容：【实验方案设计】我们分别运用Euler法，改进Euler法，RK方法和Adams隐式方法对同一问题进行求解，将数值解和解析解画在同一图像中，比较数值解的精度大小，得出结论。
【实验过程】（实验步骤、记录、数据、分析）我们首先来回顾一下原题：对于给定初值问题：1. 求出其解析解并用Matlab画出其图形；2. 采用Euler法取步长为0.5和0.25数值求解（2.16）,并将结果画在同一幅图中，比较两者精度;3. 采用改进Euler法求解（2.16），步长取为0.5；4. 采用四级Runge-Kutta法求解（2.16），步长取为0.5；5. 采用Adams四阶隐格式计算（2.16），初值可由四级Runge-Kutta格式确定。
下面，我们分五个步骤来完成这个问题：步骤一，求出（2.16）式的解析解并用Matlab 画出其图形； ，用Matlab 做出函数在上的图像，见下图：00.51 1.52 2.53 3.54 4.550.511.522.533.5x 1015y=exp(1/3 t 3-1.2t)exact solution图一 初值问题的解析解的图像步骤二，采用Euler 法取步长为0.5和0.25数值求解(2.16),并将结果画在同一幅图中，比较两者精度;我们采用Euler 法取步长为0.5和0.25数值求解，并且将数值解与解析解在一个图中呈现，见下图：00.51 1.52 2.53 3.54 4.550.511.522.533.5x 1015Numerical solution of Euler and exact solutionexact solution h=0.5h=0.25图二 Euler 方法的计算结果与解析解的比较从图像中不难看出，采用Euler 法取步长为0.5和0.25数值求解的误差不尽相同，也就是两种方法的计算精度不同，不妨将两者的绝对误差作图，可以使两种方法的精度更加直观化，见下图：00.51 1.52 2.53 3.54 4.550.511.522.533.5x 1015Absolute error of numerical solution and exact solutionh=0.5h=0.25图三 不同步长的Euler 法的计算结果与解析解的绝对误差的比较 从图像中我们不难看出，步长为0.25的Euler 法比步长为0.5的Euler 法的精度更高。
常微分方程数值解法及其应用的开题报告
常微分方程数值解法及其应用的开题报告开题报告题目：常微分方程数值解法及其应用一、研究背景及意义常微分方程（ODE）在各种自然科学和工程学科中都有着广泛应用，例如力学、物理学、化学、经济学、生物学等。
由于许多ODE无法精确求解，因此需要求解数值解。
数值解法的研究是计算数学的一个重要方向，它是许多实际问题求解的基础。
在实际应用中，我们经常遇到大型ODE系统的数值求解。
因此，如何设计高效的数值方法来求解ODE是非常重要的研究方向。
另外，ODE数值解法和应用也成为了计算数学、计算机科学、物理学等领域交叉、融合的重要部分。
二、研究现状及问题随着计算机技术的飞速发展，ODE数值解法也在不断更新和发展。
常用的ODE数值解法主要有欧拉法、龙格-库塔法、多步法等。
这些方法具有一定的理论基础和实用价值，但仍存在不足之处。
一些ODE模型的求解仍然是一项具有挑战性的任务，例如具有强非线性和奇点特性的方程。
此外，针对大型ODE系统求解的算法和并行计算方法也需要进一步研究和优化。
三、研究内容和方法本文将从ODE数值解法和应用的角度来进行研究，主要包括以下内容：1. ODE数值解法：对欧拉法、龙格-库塔法、多步法等常见数值方法进行深入研究和比较，并探讨其稳定性、精度和复杂度等方面的特点。
2. 非线性ODE求解：对具有非线性特点的ODE进行求解，探讨如何利用牛顿法、拟牛顿法等优化算法加速求解过程；3. 大型ODE系统求解：探讨如何利用稀疏矩阵技术和并行计算技术加速大型ODE系统的求解；4. 应用实例：在物理、生物领域等应用实例中，探讨ODE数值解法的应用和实现过程。
本文所采用的研究方法包括理论分析、数值实验和仿真模拟等。
四、预期成果和贡献1. 系统阐述ODE数值解法，为不同应用领域提供可供选择的计算方法和应用实例。
2. 针对具有强非线性和奇点特性的方程，提出了利用数值优化算法求解ODE的方法，实现更高效稳定的求解。
3. 探讨稀疏矩阵技术和并行计算技术在大型ODE系统求解中的应用，实现更快速、高效的求解。
常微分方程数值解实验报告
常微分方程数值解实验报告实验报告：常微分方程数值解1.引言常微分方程（Ordinary Differential Equations, ODEs）是数学领域中一个重要的研究对象，涉及到许多自然科学和工程技术领域的问题。
解常微分方程的数值方法是一种求解差分方程的方法，通过计算机找到方程的近似解，对于模拟和预测连续过程非常有用。
本实验旨在通过数值解法，验证和应用常微分方程的解，并比较不同数值方法的精度和效率。
2.实验目的2.1理解常微分方程的基本概念和数值解法；2.2掌握将常微分方程转化为数值求解问题的基本方法；2.3运用数值解法求解常微分方程；2.4比较不同数值解法的精度和效率。
3.实验内容3.1 欧拉方法（Euler Method）给定一个一阶常微分方程dy/dx=f(x,y)，通过将其离散为差分形式，欧拉方法可以通过以下递推公式来求解：y_{n+1}=y_n+h*f(x_n,y_n)其中，h为步长，x_n和y_n为当前的x和y值。
3.2 改进的欧拉方法（Improved Euler Method）改进的欧拉方法使用欧拉方法的斜率的平均值来估计每一步中的斜率。
具体公式如下：k1=f(x_n,y_n)k2=f(x_n+h,y_n+h*k1)y_{n+1}=y_n+h*((k1+k2)/2)3.3 二阶龙格-库塔法（Second-order Runge-Kutta Method）二阶龙格-库塔法通过计算每个步骤中的两个斜率来估计每个步长中的斜率。
具体公式如下：k1=f(x_n,y_n)k2=f(x_n+h/2,y_n+(h/2)*k1)y_{n+1}=y_n+h*k24.实验步骤4.1选取常微分方程，并将其转化为数值求解问题的形式；4.2根据给定的初始条件和步长，使用欧拉方法、改进的欧拉方法和二阶龙格-库塔法求解该方程；4.3比较三种方法的数值解与理论解的差异，并分析其精度和效率；4.4尝试不同的步长，观察相应的数值解的变化。
微分方程数值解实验报告
微分方程数值解实验报告实验目的：掌握微分方程数值解的基本方法，能够利用计算机编程求解微分方程。
实验原理：微分方程是自然科学与工程技术中常见的数学模型，它描述了变量之间的关系及其随时间、空间的变化规律。
解微分方程是研究和应用微分方程的基础，但有很多微分方程无法找到解析解，只能通过数值方法进行求解。
本实验采用欧拉方法和改进的欧拉方法求解微分方程的初值问题：$$\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}误差-步长曲线如下：实验结论：通过对比欧拉方法和改进的欧拉方法的数值解与解析解的误差，可以发现改进的欧拉方法具有更高的精度和收敛性。
实验报告七常微分方程初值问题的数值解法
浙江大学城市学院实验报告课程名称数值计算方法实验项目名称常微分方程初值问题的数值解法 实验成绩指导老师签名日期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四.实验结果与分析。
实验5常微分方程的数值解
实验5 常微分方程的数值解概要：将装满放射性废物的圆桶扔到水深300ft 的海底，圆桶体积55gal ，装满废料的桶重为527.436lbf ，在海中浮力为470.327lbf 。
此外，下沉时受到的阻力与速度成正比，比例系数为0.08lbf/s 。
实验发现当圆桶速度超过40ft/s 时，就会因与海底冲撞而破裂。
要求：（1）建立解决上述问题的微分方程模型（2）用数值和解析两种方法求解微分方程，并回答谁赢得了官司。
模型建立由牛顿第二定律可列出圆桶下沉速度的微分方程：dv G F f G F bv dt m m ----==其中G 为圆桶重量，F 为浮力，b 为下沉阻力与速度关系的比例系数。
换算到国际单位制，dept=300*0.3048=91.4400 海深（m ）ve=40*0.3048=12.1920 速度极限（超过就会使圆筒碰撞破裂）(m/s) G=527.436*0.4536*9.8＝2344.6 圆筒重量（N ） F=470.327*0.4536*9.8＝2090.7 浮力(N)m=527.436*0.4536＝239.24 圆筒质量（kg ） b=0.08*0.4536*9.8/0.3048＝1.1667 比例系数(Ns/m) 模型求解 一．求数值解Matlab 程序如下： sd.m:function dx=sd(t,x,G,F,m,b) dx=[(G-F-b*x)/m];%微分方程sddraw.m: clear;G=527.436*0.4536*9.8;%圆筒重量（N ） F=470.327*0.4536*9.8;%浮力(N)m=527.436*0.4536;%圆筒质量（kg ）b=0.08*0.4536*9.8/0.3048%比例系数(Ns/m) h=0.1;%所取时间点间隔ts=[0:h:2000];%粗略估计到时间2000 x0=0;%初始条件opt=odeset('reltol',1e-3,'abstol',1e-6);%相对误差1e-6，绝对误差1e-9 [t,x]=ode45(@sd,ts,x0,opt,G ,F,m,b);%使用5级4阶龙格—库塔公式计算 %[t,x]%输出t,x(t),y(t)plot(t,x,'-'),grid%输出v(t)的图形 xlabel('t'); ylabel('v(t)');%用辛普森公式对速度积分求出下沉深度 T=20;%估计20s 以内降到海底for i=0:2:10*T%作图时间间隔为0.2 y=x(1:(i+1)); k=length(y);a1=[y(2:2:k-1)];s1=sum(a1); a2=[y(3:2:k-1)];s2=sum(a2);z4((i+2)/2)=(y(1)+y(k)+4*s1+2*s2)*h/3;%辛普森公式求深度 endi=[0:2:10*T]; figure;de=300.*0.3048.*ones(5*T+1,1);%海深ve=40.*0.3048*[1 1];%速度极限值（超过就会使圆筒碰撞破裂）plot(x(i+1),z4',x(i+1),de,ve,[0 z4(5*T+1)]);%作出速度－深度图线，同时画出海深和速度要求grid;gtext('dept'),gtext('Vmax');xlabel('v');ylabel('dept(v)');figure;plot(i/10,z4');%作出时间－下降深度曲线grid;xlabel('t');ylabel('dept(t)');求解结果如下图：速度—时间曲线：可以看到经过足够长的时间后，如若桶没有落到海底，它的速度会趋于常值。
微分方程数值解法实验报告
微分方程数值解法实验报告姓名： 班级： 学号：一：问题描述求解边值问题：()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、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
常微分方程数值解实验报告学院:数学与信息科学专业:信息与计算科学姓名：郑思义学号:2０121６52４课程：常微分方程数值解实验一：常微分方程得数值解法1、分别用Eｕleｒ法、改进得Euler法(预报校正格式)与S—K法求解初值问题。
(h＝0、1)并与真解作比较。
1、1实验代码:%欧拉法functioｎ [ｘ，y]=naeuler(dｙｆｕｎ，xspan,y０,ｈ)%ｄyfun就是常微分方程,xｓpan就是x得取值范围，ｙ０就是初值,h就是步长x＝ｘsｐan(1）:h:xspan(2)；y(1）=y０；ｆｏrｎ=1:length(x)-１ｙ（n+１）=ｙ(n)+h*feval（dｙｆun,ｘ(ｎ),ｙ(ｎ））;eｎd％改进得欧拉法ｆunction [x,ｍ,ｙ]=naｅuler2(ｄyfun,xspaｎ,y０,h)%ｄyfｕn就是常微分方程,xｓpan就是ｘ得取值范围,y0就是初值，h就是步长。
%返回值ｘ为x取值,ｍ为预报解,ｙ为校正解x＝xspan(1)：h:xspaｎ(2);y（1)=y0;m=zerｏs（lenｇth(x)－1，１);foｒｎ=1:leｎgｔh(x）-1k1＝fｅｖal（dｙfun,x(n),y(n）)；y(n＋1)=y（n)+h*k1;ｍ（n)=y（n+1);k2＝ｆevａｌ(dyfuｎ,x(n＋1）,y(n+1）);y（ｎ+1）=y(n)+h*（k１+k2)／2;end%四阶S—K法functioｎ [ｘ，ｙ]=rk(ｄyfun,xspａn,y0,h）%dｙｆun就是常微分方程,xｓpaｎ就是ｘ得取值范围,y0就是初值，h就是步长。
ｘ＝xsｐan(1):h：xsｐan(2)；ｙ(１)=y0;for n=1:lｅｎgth(x)-１ｋ1=fｅｖａl(dyfｕn，x(n），y(ｎ）)；k2=fｅｖal(dyfun,x(ｎ)+ｈ/2，y(ｎ）+(h＊ｋ1)/2);k3=feval(ｄyfuｎ,ｘ(n)+ｈ/２，y(n)+(h*ｋ2)/2);k4=ｆevａl(dｙfun，x(ｎ)+h,y(n)+ｈ＊k3）；y(n+１)＝ｙ(n)+（h/6）*(k1+２*k2+2*k3+k４）;enｄ％主程序x=［０:0、1:１］;y＝ｅxp(-ｘ)+x；dyfun=inline('-ｙ+x+1'）；[x1，y1]=naeuler(ｄyｆｕn,[0,１]，1,0、1);［x2,m，y2]＝naeuler2（dｙfuｎ,[０，1］,1,０、1）;[x3,y3]=rｋ(dyｆuｎ,[0,１］,1,０、１);ploｔ(x，y,'r',x1,y１,'+'，x2,y2,＇*',x3，y３,'o');xlabｅl('ｘ');yｌabel('y＇);legｅnd(＇y为真解',＇y1为欧拉解＇,'ｙ2为改进欧拉解','y3为S—K解','Locatiｏn','ＮorthWeｓｔ＇)；1、２实验结果：x 真解y 欧拉解y1 预报值ｍ校正值y2 S—K解y30、０１、0000 1、0000 1、０000 １、000０0、1 1、0048 1、000０1、００00 １、005０１、0０4８0、2 １、0187 １、0100 1、０145 1、０190 1、0187０、3 1、0408 １、02９0 1、０３7１1、0412 1、0４0８０、4 １、0７03 1、0５６1 1、06７1 1、0708 １、0７030、5 1、10６5 1、09０５1、1０3７1、107１1、１0６50、6 1、１488 1、１３1４１、1464 １、149４1、１4880、7 1、1966 １、17８3 1、194５1、19７2 1、１9660、８1、2493 １、2３05 1、2475 1、2500 1、２４９30、9 1、306６1、２87４1、３050 1、3072 １、３0661、0 1、3679 1、34８７1、36６5 1、３６85 1、36７92、选取一种理论上收敛但就是不稳定得算法对问题1进行计算,并与真解作比较。
(选改进得欧拉法）２、１实验思路:算法得稳定性就是与步长h密切相关得。
而对于问题一而言，取定步长h=0、1不论就是单步法或低阶多步法都就是稳定得算法。
所以考虑改变h取值范围,借此分析不同步长会对结果造成什么影响。
故依次采用h=２、0、2、２、２、4、2、６得改进欧拉法。
2、2实验代码:%%主程序x=[0：３:30］;y＝exp(－ｘ)+ｘ;dyｆuｎ=inliｎｅ('-y＋ｘ+１');[x1,ｍ1,y１］=naeｕler2（ｄｙfｕn,［０,20],1，2);[x2，m２，y2]=naeｕｌeｒ２（dｙfｕｎ，[0,２2],１，2、2);[x3，ｍ3,y3］＝naeuｌｅr2(dｙfun,[0,２4]，1,２、４);[x4,m４,ｙ4］=naｅｕlｅr2（dyfun,[0,2６],1,２、6)；subｐｌｏt(2，2,１）pｌｏｔ(x,y,'r'，x1,y１，＇+');ｘlabel('h=２、0')；subplｏt（2,2,２)plot(x,ｙ，＇r＇,ｘ2，y2，'+'）；xlａbel('h=２、２');sｕｂpｌot(2,２，３)ploｔ(x,ｙ,＇ｒ'，x3,ｙ3,'+＇)；xlabeｌ(＇h=2、4＇）;sｕbplot(２，2,4)pｌot(ｘ,y,'r＇,x4，y４,'＋');xlaｂel('h=２、6');2、3实验结果:x h=２、0 h=2、2 h=２、4 h=2、60、０1、0０00 １、0000 1、0000 １、０0000、１3、０000 3、4200 3、８８00 ４、38０00、2 5、000０5、8884 6、９90４8、368４0、3 7、0０00 ８、４1５8 10、44１8 1３、4398０、4 9、000０11、0１53 １4、3979 20、43880、5 11、0000 1３、70２7 １９、1０08 30、８6900、6 13、0000 16、4973 ２４、909２４７、4０68０、７15、０0０0 １9、4227 32、３5３6 74、８１610、8 17、０0０0 22、５０７７42、2194 1２1、57670、９19、000０25、7874 55、66８7 202、7８251、0 ２１、0００0 29、30４6 ７4、4２１7 345、３０0８实验结果分析:从实验1结果可以瞧出,在算法满足收敛性与稳定性得前提下,Elｕeｒ法虽然计算并不复杂,凡就是精度不足,反观改进得Ｅlｕer法与S—K法虽然计算略微复杂但就是结果很精确。
实验２改变了步长,导致算法理论上收敛但就是不满足稳定性。
结果表示步长ｈ越大,结果越失真。
对于同一个问题,步长h得选取变得尤为重要,这三种单步法得绝对稳定区间并不一样,所以并没有一种方法就是万能得,我们应该根据不同得步长来选取合适得方法。
实验二:Ｒｉtz-Gａlerｋｉｎ方法与有限差分法1、用中心差分格式求解边值问题取步长ｈ=０、1，并与真解作比较。
1、1实验代码:%中心差分法fｕｎction U=fdｍ(xspan,y0，y１,ｈ)%xsｐａn为x取值范围,y0,y1为边界条件,h为步长N=１/ｈ;d=zeroｓ（1,N－1）;foｒ i=1:Nx(i)=ｘspan(1)＋i＊ｈ;q(i)=１;ｆ（i）＝x（i)；endfoｒｉ=1:Ｎ-１d(i)=q(i)＊h*h＋２；eｎda=diaｇ(d)；b=zｅros（N-1）;ｃ＝zeｒos(N-１);ｆｏr i=１：N-2b(i+1,i)=-1;enｄfoｒi＝1:Ｎ-2ｃ(i,i+1）=－1;endA=ａ+ｂ+c；ｆoｒｉ=2：N－2B(i,1)=ｆ（i)*h＊ｈ;eｎdＢ（1,１)=ｆ(１)＊h＊ｈ＋y０;B(N-1,1)＝f(N-1)＊ｈ*h+y1;Ｕ= iｎv(A)*B;%主程序x=0:0、1:1;ｙ＝ｘ+(exp(1)*exp(－x))/(exｐ(２)－1)-(exp(1)*eｘｐ(x）)／(ｅｘｐ(2)－1); y1=ｆdm([０,1］,0,0，0、１);ｙ1=［０,y1',0];plot（ｘ,y,'r'，x,y1，'+＇）ｘlabｅl('x');ylabｅl（'y＇)；leｇend（＇y真解'，'y1中心差分法＇,'Location',＇NｏｒｔhＷest＇)；1、2实验结果:x ｙ真解y1中心差分法0、0 0、0000 ０、0００0０、1 0、０１4８０、0148０、２0、028７0、02８70、3 0、04０９0、０４0８0、4 0、０5０5 0、05040、5 0、０56６0、05650、6 ０、0583 ０、０５820、7 0、0545 0、05４5０、8 0、０4４3 0、０44３0、9 0、０265 ０、０265１、０0、0000 0、00002、用Rｉtz-Galeｒｋｉｎ方法求解上述问题,并与真值作比较,列表画图。
２、１实验代码:%Rｉｔz_Galｅrｋｉｎ法funcｔion vu=Ｒitz_Gａlerkin(x0,y0，x1，y１,h)%x0,ｘ1为x取值范围，y0，y1为边界条件，h为步长N=1／h;ｓyms x;ｆor i=1:Ｎfａi（ｉ）=x*(１-ｘ)＊(x^(i-1));dfai(i)=diff(x*（1-ｘ）＊（x^(i－1))）;ｅnｄforｉ=1：Nfor j=1：Nfun=ｄfai(i）*dｆaｉ（ｊ)+fai(ｉ)*fai(j)；A（ｉ,ｊ）=int(fuｎ,x,0,1）;ｅndfｕn=x*faｉ(i)+dfai(i）；ｆ(i)＝int(ｆｕn,x,0，1);eｎdc=ｉnv(A)*f＇;product＝c、*faｉ＇；sｕm=0;fｏr i=１:Ｎsuｍ＝sｕm＋product(i）;ｅndvu=[];for y＝0:h:1v=sｕbｓ(sum，ｘ,y)；ｖu=［vｕ,ｖ];endy=0：h:１；ｙy=０：0、１:1;ｕ=sｉｎ(ｙｙ)/siｎ(1)-ｙy；u＝vｐａ(u,５);vu=ｖpａ(vu,５);%主程序ｘ＝0:0、1:1;ｙ=x+(exｐ(1)＊exｐ(-ｘ）)/(exｐ(２)-1)-(exｐ(1)＊exp(x）)/(ｅｘp(2)－1); y1＝Rｉtｚ_Galerkin(0,0,１,0,０、1);y１=ｄoublｅ(y１）;ｐloｔ(ｘ,y,'ｒ',x，y1,'+＇）xlaｂel(＇x');ylabｅl('y');ｌegeｎｄ（'ｙ为真解',＇ｙ1为Ｒ—G法'，'Loｃation','NorthWest')；２、2实验结果:x y真解ｙ１Ｒ—Ｇ法０、0 ０、０000 0、0０000、１０、０１48 ０、０1480、２0、０28７0、0２870、3 0、０409 0、04090、4 0、0５05 ０、０5050、5 0、０566 0、05660、6 0、0５83 0、０５8３0、7 0、0５4５0、054５0、８０、0443 0、0４43０、9 0、０265 0、0２６51、0 ０、00０0 ０、00003、若边值条件为y(0)=0,ｙ（１)=1;则上述问题得数值解法怎样变化？结果如何?程序运算出来真解与数值解完全一样。