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

合集下载

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

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

常微分方程数值解实验报告学院:数学与信息科学专业:信息与计算科学姓名:郑思义学号: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进行计算,并与真解作比较。

数值分析第九章常微分方程数值解法

数值分析第九章常微分方程数值解法
高斯-赛德尔迭代法
松弛法
通过迭代更新函数值并逐步放松约束 条件来逼近解,适用于刚性和非刚性 问题。
利用线性组合迭代函数值来逼近解, 具有更高的收敛速度和稳定性。
03
数值解法的稳定性分析
数值解法的稳定性定义
数值解法的稳定性是指当微分方程的初值有微小的扰动时, 其数值解的近似值的变化情况。如果数值解在微小扰动下变 化较小,则称该数值方法是稳定的。
更高的精度和稳定性。
数值逼近法
泰勒级数法
将微分方程的解展开为泰勒级数,通过截断级数来逼 近解。
多项式逼近法
利用多项式来逼近微分方程的解,通过选取合适的基 函数和系数来提高逼近精度。
样条插值法
利用样条函数来逼近微分方程的解,具有更好的光滑 性和连续性。
迭代法
雅可比迭代法
通过迭代更新函数值来逼近微分方程 的解,具有简单易行的优点。
初值和边界条件的处理
根据实际问题,合理设定初值和边界 条件,以获得更准确的数值解。
收敛性和误差分析
对数值解进行收敛性和误差分析,评 估解的精度和稳定性。
数值解法的应用案例分析
人口增长模型
通过数值解法求解人口增长模型,预测未来人口数量,为政策制 定提供依据。
化学反应动力学
利用数值解法研究化学反应的动力学过程,模拟反应过程和结果。
数值分析第九章常微分方 程数值解法
• 引言 • 常微分方程数值解法的基本思想 • 数值解法的稳定性分析 • 数值解法的收敛性和误差分析 • 数值解法的实现和应用案例
01
引言
常微分方程的应用背景
自然科学
描述物理、化学、生物等自然 现象的变化规律。
工程领域
控制系统设计、航天器轨道计 算等。

实验报告——常微分方程的数值解法

实验报告——常微分方程的数值解法

实验报告实验项目名称常微分方程的数值解法实验室数学实验室所属课程名称微分方程数值解实验类型上机实验实验日期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。

数值分析常微分方程数值解法

数值分析常微分方程数值解法
7
第8页/共105页
➢ 数值积分方法(Euler公式)
设将方程 y=f (x, y)的两端从 xn 到xn+1 求积分, 得
y( xn1) y( xn )
xn1 f ( x, y( x))dx :
xn
xn1 F ( x)dx
xn
用不同的数值积分方法近似上式右端积分, 可以得到计算 y(xn+1)的不同的差分格 式.
h2 2
y''( )
Rn1
:
y( xn1)
yn1
h2 2
y''( )
h2 2
y''( xn ) O(h3 ).
局部截断误差主项
19
第20页/共105页
➢ 向后Euler法的局部截断误差
向后Euler法的计算公式
yn1 yn hf ( xn1, yn1 ), n 0, 1, 2,
定义其局部截断误差为
y 计算 的n递1 推公式,此类计算格式统称为差分格式.
3
第4页/共105页
数值求解一阶常微分方程初值问题
y' f ( x, y), a x b,
y(a)
y0
难点: 如何离散 y ?
➢ 常见离散方法
差商近似导数 数值积分方法 Taylor展开方法
4
第5页/共105页
➢ 差商近似导数(Euler公式)
(0 x 1)
y(0) 1.
解 计算公式为
yn1
yn
hfn
yn
h( yn
2xn ), yn
y0 1.0
n 0, 1, 2,
取步长h=0.1, 计算结果见下表
13

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

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

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

常微分方程数值解实验

常微分方程数值解实验
刚性
多步法,Gear’s反向
数值积分,精度中等
若ode45失效时,
可尝试使用
ode23s
刚性
一步法,2阶Rosebrock算法,
低精度。
当精度较低时,
计算时间比ode15s短
odefx为显式常微分方程 中的 ,t为求解区间,要获得问题在其他指定点 上的解,则令t=[t0,t1,t2,…](要求 单调),y0初始条件。
MATLAB 中有几个专门用于求解常微分方程的函数,它们的设计思想基于Runge-Kutta方法,基本设计思想为:从改进的欧拉方法比欧拉方法精度高的缘由着手,如果在区间[ x1, xi+1]多取几个点的斜率值,然后求取平均值,则可以构造出精度更高的计算方法。 这些函数主要包括:ode45、ode23、ode15s、ode113、ode23s、ode23t、ode23tb. 其中最常用的是函数ode45,该函数采用变步长四阶五阶Runge-Kutta法求数值解,并采用自适应变步长的求解方法。ode23采用二阶三阶Runge-Kutta法求数值解,与ode45类似,只是精度低一些。ode15s用来求刚性方程组。
43
4月22日
588
666
28
46
4月23日
693
782
35
55
4月24日
774
863
39
64
4月25日
877
954
42
73
4月26日
988
1093
48
76
4月27日
1114
1255
56
78
4月28日
1199
1275
59
78
4月29日

求常微分方程的数值解

求常微分方程的数值解

求常微分方程的数值解一、背景介绍常微分方程(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)$。

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

数学与计算科学学院实验报告
实验项目名称常微分方程数值解
所属课程名称数值方法B
实验类型验证
实验日期 2013.11.11
班级
学号
姓名
成绩
【实验过程】(实验步骤、记录、数据、分析)
注:以下图形是通过Excel 表格处理数据得出,并未通过MATLAB 编程序所得!
1、1(0)1dy
y x dx y ⎧=-++⎪⎨⎪=⎩
由题可知精确解为:x y x e -=+,当x=0时,y(x)=0。

h=0.1
表1 h=0.1时三个方法与精确值的真值表
图1 h=0.1时三个方法走势图
步长 Euler 法 预估校正法 经典四阶库 精确值 0.1 1.010000 1.005000 1.004838 1.249080 0.2 1.029000 1.019025 1.018731 1.055455 0.3 1.056100 1.041218 1.040818 1.091217 0.4 1.090490 1.070802 1.070320 1.131803 0.5 1.131441 1.107076 1.106531 1.176851 0.6 1.178297 1.149404 1.148812 1.226025 0.7 1.230467 1.197211 1.196586 1.279016 0.8 1.287421 1.249975 1.249329 1.335536 0.9 1.348678 1.307228 1.306570 1.395322 1.0
1.413811
1.368541
1.367880
1.458127
h=0.05(此时将源程序中i的围进行扩大,即for(i=0;i<20;i++))
表2 h=0.05时三个方法与精确值的真值表步长Euler法预估校正法经典四阶库精确值
0.05 1.002500 1.001250 1.001229 1.011721
0.10 1.007375 1.004877 1.004837 1.024908
0.15 1.014506 1.010764 1.010708 1.039504
0.20 1.023781 1.018802 1.018731 1.055455
0.25 1.035092 1.028885 1.028801 1.072710
0.30 1.048337 1.040915 1.040818 1.091217
0.35 1.063421 1.054795 1.054688 1.110931
0.40 1.080250 1.070436 1.070320 1.131801
0.45 1.098737 1.087752 1.087628 1.153791
0.50 1.118800 1.106662 1.106531 1.176851
0.55 1.140360 1.127087 1.126950 1.200942
0.60 1.163342 1.148954 1.148812 1.226025
0.65 1.187675 1.172193 1.172046 1.252062
0.70 1.213291 1.196736 1.196585 1.279016
0.75 1.240127 1.222520 1.222367 1.306852
0.80 1.268121 1.249485 1.249329 1.335536
0.85 1.297215 1.277572 1.277415 1.365037
0.90 1.327354 1.306728 1.306570 1.395322
0.95 1.358486 1.336900 1.336741 1.426362
1.00 1.390562 1.368039 1.367880 1.458127
图2 h=0.05时三个方法走势图
2、(0)1x
dy xe y dx y -⎧=-⎪⎨⎪=⎩
由题可知精确解为:21
(2)2
x x e -+,当x=0时,y(x)=0。

h=0.1
表3 h=0.1时三个方法与精确值的真值表
步长 Euler 法 预估校正法 经典四阶库 精确值 0.1 0.900000 0.909625 0.909428 0.929533 0.2 0.819249 0.835927 0.835593 0.872564 0.3 0.754433 0.776081 0.775655 0.826822 0.4 0.702726 0.727671 0.727189 0.790348 0.5 0.661726 0.688636 0.688127 0.761457 0.6 0.629396 0.657225 0.656711 0.738709 0.7 0.604018 0.631957 0.631453 0.720874 0.8 0.584147 0.611582 0.611100 0.706908 0.9 0.568575 0.595050 0.594599 0.695927 1.0 0.556297 0.581487 0.581072
0.687191
图3 h=0.1时三个方法走势图
h=0.05(此时将源程序中i的围进行扩大,即for(i=0;i<20;i++))
表4 h=0.05时三个方法与精确值的真值表
步长Euler法预估校正法经典四阶库精确值0.050.9500000.9524520.9524270.962924 0.100.9049040.9094740.9094280.929533 0.150.8642840.8706700.8706060.899511 0.200.8277410.8356710.8355920.872564 0.250.7949080.8041370.8040470.848419 0.300.7654470.7757550.7756550.826822 0.350.7390430.7502320.7501250.807538 0.40 0.7154070.7273020.7271890.790348 0.45 0.6942720.7067150.7065990.775050 0.50 0.6753940.6882450.6881260.761457 0.55 0.6585460.6716820.6715610.749397 0.60 0.6435190.6568300.6567100.738709 0.65 0.6301240.6435140.6433950.729247 0.70 0.6181850.6315700.6314530.720874 0.75 0.6075410.6208480.6207330.713466 0.80 0.5980460.6112110.6111000.706908 0.85 0.5895650.6025350.6024260.701094 0.90 0.5819760.5947030.5945990.695927
0.95 0.5751670.5876120.5875120.691320
1.00 0.5690350.5811670.5810710.687191
图4 h=0.05时三个方法走势图
附录1:源程序
附录2:实验报告填写说明
1.实验项目名称:要求与实验教学大纲一致。

2.实验目的:目的要明确,要抓住重点,符合实验教学大纲要求。

3.实验原理:简要说明本实验项目所涉及的理论知识。

4.实验环境:实验用的软、硬件环境。

5.实验方案(思路、步骤和方法等):这是实验报告极其重要的容。

概括整个实验过程。

对于验证性实验,要写明依据何种原理、操作方法进行实验,要写明需要经过哪几个步骤来实现其操作。

对于设计性和综合性实验,在上述容基础上还应该画出流程图、设计思路和设计方法,再配以相应的文字说明。

对于创新性实验,还应注明其创新点、特色。

6.实验过程(实验中涉及的记录、数据、分析):写明具体实验方案的具体实施步骤,包括实验过程中的记录、数据和相应的分析。

7.实验结论(结果):根据实验过程中得到的结果,做出结论。

8.实验小结:本次实验心得体会、思考和建议。

9.指导教师评语及成绩:指导教师依据学生的实际报告容,给出本次实验报告的评价。

相关文档
最新文档