数值分析实验报告之常微分方程数值解

常微分方程数值解实验报告学院：数学与信息科学专业：信息与计算科学姓名：郑思义学号：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 法的精度更高。

许多实际问题的数学模型是微分方程或微分方程的定解问题。

如物体运动、电路振荡、化学反映及生物群体的变化等。

常微分方程可分为线性、非线性、高阶方程与方程组等类；线性方程包含于非线性类中，高阶方程可化为一阶方程组。

若方程组中的所有未知量视作一个向量，则方程组可写成向量形式的单个方程。

因此研究一阶微分方程的初值问题⎪⎩⎪⎨⎧=≤≤=0)(),(y a y bx a y x f dxdy， (9-1) 的数值解法具有典型性。

常微分方程的解能用初等函数、特殊函数或它们的级数与积分表达的很少。

用解析方法只能求出线性常系数等特殊类型的方程的解。

对非线性方程来说，解析方法一般是无能为力的，即使某些解具有解析表达式，这个表达式也可能非常复杂而不便计算。

因此研究微分方程的数值解法是非常必要的。

只有保证问题(9-1)的解存在唯一的前提下，研究其数值解法或者说寻求其数值解才有意义。

由常微分方程的理论知，如果(9-1)中的),(y x f 满足条件（1）),(y x f 在区域} ),({+∞<<∞-≤≤=y b x a y x D ，上连续； （2）),(y x f 在上关于满足Lipschitz 条件，即存在常数，使得y y L y x f y x f -≤-),(),(则初值问题(9-1)在区间],[b a 上存在惟一的连续解)(x y y =。

在下面的讨论中，我们总假定方程满足以上两个条件。

所谓数值解法，就是求问题(9-1)的解)(x y y =在若干点b x x x x a N =<<<<= 210处的近似值),,2,1(N n y n =的方法。

),,2,1(N n y n =称为问题(9-1)的数值解，n n x x h -=+1称为由到1+n x 的步长。

今后如无特别说明，我们总假定步长为常量。

建立数值解法，首先要将微分方程离散化，一般采用以下几种方法： (1) 用差商近似导数在问题(9-1)中，若用向前差商hx y x y n n )()(1-+代替)(n x y '，则得)1,,1,0( ))(,()()(1-=≈-+N n x y x f hx y x y n n n n n)(n x y 用其近似值代替，所得结果作为)(1+n x y 的近似值，记为1+n y ，则有 1(,) (0,1,,1)n n n n y y hf x y n N +=+=-这样，问题(9-1)的近似解可通过求解下述问题100(,) (0,1,,1)()n n n n y y hf x y n N y y x +=+=-⎧⎨=⎩(9-2)得到，按式(9-2)由初值经过步迭代，可逐次算出N y y y ,,21。

常微分方程数值解实验报告实验报告：常微分方程数值解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尝试不同的步长，观察相应的数值解的变化。

求常微分方程的数值解一、背景介绍常微分方程（Ordinary Differential Equation，ODE）是描述自然界中变化的数学模型。

常微分方程的解析解往往难以求得，因此需要寻找数值解来近似地描述其行为。

求解常微分方程的数值方法主要有欧拉法、改进欧拉法、龙格-库塔法等。

二、数值方法1. 欧拉法欧拉法是最简单的求解常微分方程的数值方法之一。

它基于导数的定义，将微分方程转化为差分方程，通过迭代计算得到近似解。

欧拉法的公式如下：$$y_{n+1}=y_n+f(t_n,y_n)\Delta t$$其中，$y_n$表示第$n$个时间步长处的函数值，$f(t_n,y_n)$表示在$(t_n,y_n)$处的导数，$\Delta t$表示时间步长。

欧拉法具有易于实现和理解的优点，但精度较低。

2. 改进欧拉法（Heun方法）改进欧拉法又称Heun方法或两步龙格-库塔方法，是对欧拉法进行了精度上提升后得到的一种方法。

它利用两个斜率来近似函数值，并通过加权平均来计算下一个时间步长处的函数值。

改进欧拉法的公式如下：$$k_1=f(t_n,y_n)$$$$k_2=f(t_n+\Delta t,y_n+k_1\Delta t)$$$$y_{n+1}=y_n+\frac{1}{2}(k_1+k_2)\Delta t$$改进欧拉法比欧拉法精度更高，但计算量也更大。

3. 龙格-库塔法（RK4方法）龙格-库塔法是求解常微分方程中最常用的数值方法之一。

它通过计算多个斜率来近似函数值，并通过加权平均来计算下一个时间步长处的函数值。

RK4方法是龙格-库塔法中最常用的一种方法，其公式如下：$$k_1=f(t_n,y_n)$$$$k_2=f(t_n+\frac{\Delta t}{2},y_n+\frac{k_1\Delta t}{2})$$ $$k_3=f(t_n+\frac{\Delta t}{2},y_n+\frac{k_2\Delta t}{2})$$ $$k_4=f(t_n+\Delta t,y_n+k_3\Delta t)$$$$y_{n+1}=y_n+\frac{1}{6}(k_1+2k_2+2k_3+k_4)\Delta t$$三、数值求解步骤对于给定的常微分方程，可以通过以下步骤求解其数值解：1. 确定初值条件：确定$t=0$时刻的函数值$y(0)$。

微分方程数值解法实验报告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]。