文档库 最新最全的文档下载
当前位置:文档库 › 数值分析报告常微分方程的数值解法

数值分析报告常微分方程的数值解法

数值分析报告常微分方程的数值解法
数值分析报告常微分方程的数值解法

常微分方程的数值解法

一、 实验目的:

目的与要求:通过实验,熟悉常微分方程的数值解法的基本原理。掌握向前欧拉法、向后欧拉法、梯形法、改进欧拉法及三阶、四阶龙格-库塔法等基本算法。

二、 实验内容:

在下列方法中:向前欧拉法、向后欧拉法、梯形法、改进欧拉法及三阶、四阶龙格-库塔法 选择不同的三种算法求下面方程的数值解,并与准确解比较。

()2

,0100y x x y x y '?=+-?<

三、 实验要求:

1.编程实现。2.可以自拟实验题,要求同上。

四、 实验程序及结果:

1. 求原函数

DSolve[{y'[x]==x^2+x-y[x],y[0]==0},y[x],x] 求得:

2. 向前欧拉法

x[0]=0;y[0]=0;h=0.1; f[x_]:=;

x[n_]:=n*h;

f[u_,v_]:=u 2+u-v;

K1[n_]:=f[x[n-1],y[n-1]];

y[n_]:=y[n-1]+h*K1[n];

Table[{n,x[n],y[n],f[x[n]],f[x[n]]-y[n]},{n,0,10}]//N;

Print["\nn,x[n],y[n],f[x[n]],f[x[n]]-y[n]\n","

",MatrixForm[Table[{n,x[n],y[n],f[x[n]],f[x[n]]-y[n]},{n,0,10}]//N]]

Data1=Table[{x[n],f[x[n]]},{n,0,10}]//N;

Data2=Table[{x[n],y[n]},{n,0,10}]//N;

A0=Plot[f[x],{x,0,1},PlotStyle →{Red,Thick}];

A1=Graphics[{PointSize[0.02],Red,Point[Data1]}];

A2=ListLinePlot[Data2,PlotStyle →{Red,Dashed}];

A3=Graphics[{PointSize[0.02],Blue,Point[Data2]}];

Show[A0,A1,A2,A3,PlotRange →All]

3.改进欧拉法

f[x_,y_]:=x2+x-y;

f[x_]:=;

xy:={{0,1}};

h=0.1;

Do[xn=xy[[n]][[1]];yn=xy[[n]][[2]];

k1=f[xn,yn];k2=f[xn+h,yn+h*k1];

d=(k1+k2)*h/2 ;

Print[n," ","x",n,"=",xn+h," ","y",n,"=",yn+d," ","f(x",n,")=",f[xn+h],"

","f(x",n,")-y(x",n,")=",f[xn+h]-(yn+d)];

xy=Append[xy,{xn+h,yn+d}],{n,1,10}];

A1=Graphics[{PointSize[0.02],Red,Point[xy]}];

A2=ListLinePlot[xy,PlotStyle→Dashed];

A3=Plot[f[x],{x,0,1},PlotStyle→{Red,Thick}]; Show[A3,A2,A1]

4.四阶龙格-库塔法

f[x_,y_]:=x2+x-y;

f[x_]:=;;

y0=1;a=0;b=1;n=10;

h=(b-a)/n;

xx=Table[a+(i-1)h,{i,1,n+1}]//N;

y=Table[0,{i,1,n+1}];

y[[1]]=y0//N;

For[i=2,i n+1,i++,

a1=h f[xx[[i-1]],y[[i-1]]];

a2=h f[xx[[i-1]]+h/2,y[[i-1]]+a1/2];

a3=h f[xx[[i-1]]+h/2,y[[i-1]]+a2/2];

a4=h f[xx[[i-1]]+h,y[[i-1]]+a3];

y[[i]]=y[[i-1]]+1/6 (a1+2a2+2a3+a4);

Print[i-1," ","x",i-1,"=",xx[[i]]," ","y",i-1,"=",y[[i]]," ","f(x",i-1,")=",f[xx[[i]]],"

","f(x",n,")-y(x",n,")=",f[xx[[i]]]-y[[i]]];]

data=Table[{xx[[i]],y[[i]]},{i,1,n+1}];

A1=Graphics[{PointSize[0.02],Blue,Point[data]}];

A2=ListLinePlot[data];

A3=Plot[f[x],{x,0,1},PlotStyle {Red,Thick}];

Show[A3,A2,A1]

五、实验总结:

本次实验是为了熟悉常微分方程的数值解法的基本原理。掌握向前欧拉法、向后欧拉法、梯形法、改进欧拉法及三阶、四阶龙格-库塔法等基本算法。

对同一个方程求数值解,我采用了向前欧拉法、改进欧拉法和四阶龙格-库塔法,了解了不同算法的编程难度和运算效率。本次实验中又犯了粗心的错误,少打了一个“]”,使程序不能正确运行,花了不少时间查错。

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

《计算机数学基础》数值部分第五单元辅导 14 常微分方程的数值解法 一.重点内容 1. 欧拉公式: )心知1)a 儿+1 =儿 + hfg ,儿) m 1、 伙=0丄2,…川一 1) I 无=x Q +kh 局部截断误差是0(*)。 2. 改进欧拉公式: 预报一校正公式: 预报值 _v*+1 =儿+ hf (x k ,儿) - h - 校正值 y M = y k +-[f (x kt y k ) + /(x A+1, y M )] 即 儿+1 =儿+ £ "(忑'儿)+心+「儿+ hfg ,儿))] 或表成平均的形式: 儿=儿+ hfg ,儿) '儿=儿+"(无+】,儿) +K ) 改进欧拉法的局部截断误差是0(2) 3. 龙格一库塔法 二阶龙格一库塔法的局部截断误差是0(爪) 三阶龙格一库塔法的局部截断误差是0(护) 四阶龙格F 塔法公式:儿计=儿+ 2(匕+ 2心+ 2? + ?) 四阶龙格一库塔法的局部截断误差是0(爪)。 二实例 y' = — y — xv f2(0 < x < 0.6) 例1用欧拉法解初值问题{ ' ? -取步长/匸02计算过程保留 b (o )= 1 4位小数。 解/i=0.2. f (x )= —y —xy 2<,首先建立欧拉迭代格式 y*+i =儿+ hf g,y k ) = y k -hy k -hx k y ; =0.2 儿(4 - x k y k )(k = 0,1,2) K 2=f(x n +^h, yk+-hK\)t gg+舟人,>'n +y/?A3);

当k=0, xi=0.2 时,已知x()=0,y()=l,有 y(0?2)今i=0?2X l(4-0X 1)=0.8000 当k=\. M=0?4时,已知“=0?2」尸0?8,有 y(0?4)今2=0.2 X 0.8X(4-0.2X0.8)=0.614 4 当k=2, xs=0.6 时,已知x2=0.4,y2=0.6144,有 y(0?6)今3=0.2 X0.6144X (4-0.4 X 0.4613)=0.8000 「J, ,2 ?_ ZX 例2用欧拉预报一校正公式求解初值问题\y + v +V sinx=,取步长/?=0.2,计算 .y ⑴=1 y(0.2),y(0.4)的近似值,计算过程保留5位小数。 解步长力=0.2,此时/(x,y)=—y—fsiiu 欧拉预报一校正公式为: 预报值兀I = y k + hfg y k) - I J_ 校正值)3=儿+尹(忑,儿)+ fg,儿+1)] 有迭代格式] 预报值儿+] = y k 4-h(-y k -y; sin x k) =y k (0?8-0?2儿sin x k) < h 、—— 2 校止值y如]=儿 +尸[(一片一力sinxJ + LN+i-yl sin.v I+1)] ——?> =儿(°?9一0?1儿sin心)一0?1(儿+| +y;j sin心利) 、"M=0.別=1」)=1 时,Xj=1.2> 有 儿=yo(°?8-O?2yo sinx0) = 1 x (0.8-02x lsin 1) = 0.63171 y(1.2) ?= lx(0.9-0.1xlxsinl)-0.1(0.63171+0.631712sinl.2) = 0.71549 当 T xi=1.2, yi=0.71549 时,x2=1.4,有 y2 =儿(0.8-0?2儿sinXj) = 0.71549x(0.8-02x0.71549sinl.2) =0.47697 y(14) z y2 = 0.71549x(0.9-0.1x0.71549xsin 1.2)-0.1(0.47697+ 0.476972 sin 1.4) =0.52608 V = 8 — 3y 例3写出用四阶龙格一库塔法求解初值问题^ ‘的计算公式,取步长/匸0.2计 b(0) = 2 算y(0.4)的近似值。讣算过程保留4位小数。 解此处.心,刃=8 —3”四阶龙格一库塔法公式为 艰=儿 + % + 2? + 2勺 + ?) 1 h, y n+ y/?A3): 本例计算公式为: 0 2 呱严儿+三(32?+2?+心

常微分方程边值问题的数值解法

第8章 常微分方程边值问题的数值解法 引 言 第7章介绍了求解常微分方程初值问题的常用的数值方法;本章将介绍常微分方程的边值问题的数值方法。 只含边界条件(boundary-value condition)作为定解条件的常微分方程求解问题称为常微分方程的边值问题(boundary-value problem). 为简明起见,我们以二阶边值问题为 则边值问题(8.1.1)有唯一解。 推论 若线性边值问题 ()()()()()(),, (),()y x p x y x q x y x f x a x b y a y b αβ'''=++≤≤?? ==? (8.1.2) 满足 (1) (),()p x q x 和()f x 在[,]a b 上连续; (2) 在[,]a b 上, ()0q x >, 则边值问题(8.1.1)有唯一解。 求边值问题的近似解,有三类基本方法: (1) 差分法(difference method),也就是用差商代替微分方程及边界条件中的导数,最终化为代数方程求解; (2) 有限元法(finite element method);

(3) 把边值问题转化为初值问题,然后用求初值问题的方法求解。 差分法 8.2.1 一类特殊类型二阶线性常微分方程的边值问题的差分法 设二阶线性常微分方程的边值问题为 (8.2.1)(8.2.2) ()()()(),,(),(), y x q x y x f x a x b y a y b αβ''-=<

常微分方程初值问题的数值解法

第七章 常微分方程初值问题的数值解法 --------学习小结 一、本章学习体会 通过本章的学习,我了解了常微分方程初值问题的计算方法,对于解决那些很难求解出解析表达式的,甚至有解析表达式但是解不出具体的值的常微分方程非常有用。在这一章里求解常微分方程的基本思想是将初值问题进行离散化,然后进行迭代求解。在这里将初值问题离散化的方法有三种,分别是差商代替导数的方法、Taylor 级数法和数值积分法。常微分方程初值问题的数值解法的分类有显示方法和隐式方法,或者可以分为单步法和多步法。在这里单步法是指计算第n+1个y 的值时,只用到前一步的值,而多步法则是指计算第n+1个y 的值时,用到了前几步的值。通过对本章的学习,已经能熟练掌握如何用Taylor 级数法去求解单步法中各方法的公式和截断误差,但是对线性多步法的求解理解不怎么透切,特别是计算过程较复杂的推理。 在本章的学习过程中还遇到不少问题,比如本章知识点多,公式多,在做题时容易混淆,其次对几种R-K 公式的理解不够透彻,处理一个实际问题时,不知道选取哪一种公式,通过课本里面几种方法的计算比较得知其误差并不一样,,这个还需要自己在往后的实际应用中多多实践留意并总结。 二、本章知识梳理 常微分方程初值问题的数值解法一般概念 步长h ,取节点0,(0,1,...,)n t t nh n M =+=,且M t T ≤,则初值问题000 '(,),()y f t y t t T y t y =≤≤?? =?的数值解法的一般形式是 1(,,,...,,)0,(0,1,...,)n n n n k F t y y y h n M k ++==-

常微分方程数值解法

i.常微分方程初值问题数值解法 常微分方程初值问题的真解可以看成是从给定初始点出发的一条连续曲线。差分法是常微分方程初值问题的主要数值解法,其目的是得到若干个离散点来逼近这条解曲线。有两个基本途径。一个是用离散点上的差商近似替代微商。另一个是先对微分方程积分得到积分方程,再利用离散点作数值积分。 i.1 常微分方程差分法 考虑常微分方程初值问题:求函数()u t 满足 (,), 0du f t u t T dt =<≤ (i.1a ) 0(0)u u = (i.1b) 其中(,)f t u 是定义在区域G : 0t T ≤≤, u <∞上的连续函数,0u 和T 是给定的常数。我们假设(,)f t u 对u 满足Lipschitz 条件,即存在常数L 使得 121212(,)(,), [0,]; ,(,)f t u f t u L u u t T u u -≤-?∈∈-∞∞ (i.2) 这一条件保证了(i.1)的解是适定的,即存在,唯一,而且连续依赖于初值0u 。 通常情况下,(i.1)的精确解不可能用简单的解析表达式给出,只能求近似解。本章讨论常微分方程最常用的近似数值解法-差分方法。先来讨论最简单的Euler 法。为此,首先将求解区域[0,]T 离散化为若干个离散点: 0110N N t t t t T -=<< <<= (i.3) 其中n t hn =,0h >称为步长。 在微积分课程中我们熟知,微商(即导数)是差商的极限。反过来,差商就是微商的近似。在0t t =处,在(i.1a )中用向前差商 10()()u t u t h -代替微商du dt ,便得 10000()()(,())u t u t hf t u t ε=++ 如果忽略误差项0ε,再换个记号,用i u 代替()i u t 便得到 1000(,)u u hf t u -= 一般地,我们有 1Euler (,), 0,1, ,1n n n n u u hf t u n N +=+=-方法: (i.4) 从(i.1b) 给出的初始值0u 出发,由上式可以依次算出1,,N t t 上的差分解1,,N u u 。

常微分方程数值解

第四章常微分方程数值解 [课时安排]6学时 [教学课型]理论课 [教学目的和要求] 了解常微分方程初值问题数值解法的一些基本概念,如单步法和多步法,显式和隐式,方法的阶数,整体截断误差和局部截断误差的区别和关系等;掌握一阶常微分方程初值问题的一些常用的数值计算方法,例如欧拉(Euler)方法、改进的欧拉方法、龙贝-库塔(Runge-Kutta)方法、阿达姆斯(Adams)方法等,要注意各方法的特点及有关的理论分析;掌握构造常微分方程数值解的数值积分的构造方法和泰勒展开的构造方法的基本思想,并能具体应用它们导出一些常用的数值计算公式及评估截断误差;熟练掌握龙格-库塔(R-K)方法的基本思想,公式的推导,R-K公式中系数的确定,特别是能应用“标准四阶R-K公式”解题;掌握数值方法的收敛性和稳定性的概念,并能确定给定方法的绝对稳定性区域。[教学重点与难点] 重点:欧拉方法,改进的欧拉方法,龙贝-库塔方法。 难点:R—K方法,预估-校正公式。 [教学内容与过程] 4.1 引言 本章讨论常微分方程初值问题 (4.1.1) 的数值解法,这也是科学与工程计算经常遇到的问题,由于只有很特殊的方程能用解析方法求解,而用计算机求解常微分方程的初值问题都要采用数值方法.通常我们假定(4.1.1)中 f(x,y)对y满足Lipschitz条件,即存在常数L>0,使对,有 (4.1.2) 则初值问题(4.1.1)的解存在唯一. 假定(4.1.1)的精确解为,求它的数值解就是要在区间上的一组离散点 上求的近似.通常取 ,h称为步长,求(4.1.1)的数值解是按节点的顺序逐步 推进求得.首先,要对方程做离散逼近,求出数值解的公式,再研究公式的局部截

数值分析讲义线性方程组的解法

数值分析讲义 第三章线性方程组的解法 §3.0 引言 §3.1 雅可比(Jacobi)迭代法 §3.2 高斯-塞德尔(Gauss-Seidel)迭代法 §3.3 超松驰迭代法§3.7 三角分解法 §3.4 迭代法的收敛性§3.8 追赶法 §3.5 高斯消去法§3.9 其它应用 §3.6 高斯主元素消去法§3.10 误差分析 §3 作业讲评3 §3.11 总结

§3.0 引言 重要性:解线性代数方程组的有效方法在计算数学和科学计算中具有特殊的地位和作用.如弹性力学、电路分析、热传导和振动、以及社会科学及定量分析商业经济中的各种问题. 分类:线性方程组的解法可分为直接法和迭代法两种方法. (a) 直接法:对于给定的方程组,在没有舍入误差的假设下,能在预定的运算次数内求得精确解.最基本的直接法是Gauss消去法,重要的直接法全都受到Gauss消去法的启发.计算代价高. (b) 迭代法:基于一定的递推格式,产生逼近方程组精确解的近似序列.收敛性是其为迭代法的前提,此外,存在收敛速度与误差估计问题.简单实用,诱人.

§3.1 雅可比Jacobi 迭代法 (AX =b ) 1 基本思想: 与解f (x )=0 的不动点迭代相类似,将AX =b 改写为X =BX +f 的形式,建立雅可比方法的迭代格式:X k +1=BX (k )+f ,其中,B 称为迭代矩阵.其计算精度可控,特别适用于求解系数为大型稀疏矩阵(sparse matrices)的方程组. 2 问题: (a) 如何建立迭代格式? (b) 向量序列{X k }是否收敛以及收敛条件? 3 例题分析: 考虑解方程组??? ??=+--=-+-=--2.453.82102 .72103 21321321x x x x x x x x x (1) 其准确解为X *={1, 1.2, 1.3}. 建立与式(1)相等价的形式: ??? ??++=++=++=84.02.01.083.02.01.072 .02.01.02 13312321x x x x x x x x x (2) 据此建立迭代公式: ?????++=++=++=+++84 .02.01.083.02.01.072.02.01.0)(2)(1)1(3 )(3 )(1)1(23)(2)1(1k k k k k k k k k x x x x x x x x x (3) 取迭代初值0) 0(3 )0(2)0(1===x x x ,迭代结果如下表. JocabiMethodP31.cpp

常微分方程数值解法

第七章 常微分方程数值解法 常微分方程中只有一些典型方程能求出初等解(用初等函数表示的解),大部分的方程是求不出初等解的。另外,有些初值问题虽然有初等解,但由于形式太复杂不便于应用。因此,有必要探讨常微分方程初值问题的数值解法。本章主要介绍一阶常微分方程初值问题的欧拉法、龙格-库塔法、阿达姆斯方法,在此基础上推出一阶微分方程组与高阶方程初值问题的 数值解法;此外,还将简要介绍求解二阶常微分方程值问题的差分方法、试射法。 第一节 欧拉法 求解常微分方程初值问题 ?????==0 0)() ,(y x y y x f dx dy (1) 的数值解,就是寻求准确解)(x y 在一系列离散节点 <<<<

常微分方程的数值解

实验4 常微分方程的数值解 【实验目的】 1.掌握用MATLAB软件求微分方程初值问题数值解的方法; 2.通过实例用微分方程模型解决简化的实际问题; 3.了解欧拉方法和龙格-库塔方法的基本思想和计算公式,及稳定性等概念。 【实验内容】 题3 小型火箭初始重量为1400kg,其中包括1080kg燃料。火箭竖直向上发射时燃料燃烧率为18kg/s,由此产生32000N的推力,火箭引擎在燃料用尽时关闭。设火箭上升时空气阻力正比于速度的平方,比例系数为0.4kg/m,求引擎关闭瞬间火箭的高度、速度、加速度,及火箭到达最高点的时的高度和加速度,并画出高度、速度、加速度随时间变化的图形。 模型及其求解 火箭在上升的过程可分为两个阶段,在全过程中假设重力加速度始终保持不变,g=9.8m/s2。 在第一个过程中,火箭通过燃烧燃料产生向上的推力,同时它还受到自身重力(包括自重和该时刻剩余燃料的重量)以及与速度平方成正比的空气阻力的作用,根据牛顿第二定律,三个力的合力产生加速度,方向竖直向上。因此有如下二式: a=dv/dt=(F-mg-0.4v2)/m=(32000-0.4v2)/(1400-18t)-9.8 dh/dt=v 又知初始时刻t=0,v=0,h=0。记x(1)=h,x(2)=v,根据MATLAB 可以求出0到60秒内火箭的速度、高度、加速度随时间的变化情况。程序如下: function [ dx ] = rocket( t,x ) a=[(32000-0.4*x(2)^2)/(1400-18*t)]-9.8; dx=[x(2);a]; end ts=0:1:60;

x0=[0,0]; [t,x]=ode45(@rocket,ts,x0); h=x(:,1); v=x(:,2); a=[(32000-0.4*(v.^2))./(1400-18*t)]-9.8; [t,h,v,a]; 数据如下: t h v a 0 0 0 13.06 1.00 6.57 13.19 13.30 2.00 26.44 26.58 1 3.45 3.00 59.76 40.06 13.50 4.00 106.57 53.54 13.43 5.00 16 6.79 66.89 13.26 6.00 240.27 80.02 12.99 7.00 326.72 92.83 12.61 8.00 425.79 105.22 12.15 9.00 536.99 117.11 11.62 10.00 659.80 128.43 11.02 11.00 793.63 139.14 10.38 12.00 937.85 149.18 9.71 13.00 1091.79 158.55 9.02 14.00 1254.71 167.23 8.33 15.00 1425.93 175.22 7.65 16.00 1604.83 182.55 6.99 17.00 1790.78 189.22 6.36 18.00 1983.13 195.27 5.76 19.00 2181.24 200.75 5.21 20.00 2384.47 205.70 4.69 21.00 2592.36 210.18 4.22 22.00 2804.52 214.19 3.79 23.00 3020.56 217.79 3.41 24.00 3240.08 221.01 3.07 25.00 3462.65 223.92 2.77 26.00 3687.88 226.56 2.50 27.00 3915.58 228.97 2.27

数值分析实验2_求解线性方程组直接法

一 实验目的 1.掌握求解线性方程组的高斯消元法及列主元素法; 2. 掌握求解线性方程组的克劳特法; 3. 掌握求解线性方程组的平方根法。 二 实验内容 1.用高斯消元法求解方程组(精度要求为610-=ε): 1231231 233272212240x x x x x x x x x -+=??-+-=-??-+=? 2.用克劳特法求解上述方程组(精度要求为610-=ε)。 3. 用平方根法求解上述方程组(精度要求为610-=ε)。 4. 用列主元素法求解方程组(精度要求为610-=ε): 1231231 233432222325x x x x x x x x x -+=??-+-=??--=-? 三 实验步骤(算法)与结果 1. 程序代码(Python3.6): import numpy as np def Gauss(A,b): n=len(b) for i in range(n-1): if A[i,i]!=0: for j in range(i+1,n): m=-A[j,i]/A[i,i] A[j,i:n]=A[j,i:n]+m*A[i,i:n] b[j]=b[j]+m*b[i] for k in range(n-1,-1,-1): b[k]=(b[k]-sum(A[k,(k+1):n]*b[(k+1):n]))/A[k,k]

print(b) 运行函数: >>> A=np.array([[3,-1,2],[-1,2,-2],[2,-2,4]],dtype=np.float) >>> b=np.array([7,-1,0],dtype=np.float) >>> x=Gauss(A,b) 输出: 结果:解得原方程的解为x1=3.5,x2=-1,x3=-2.25 2 程序代码(Python3.6): import numpy as np A=np.array([[3,-1,2],[-1,2,-2],[2,-2,4]],dtype=float) L=np.array([[1,0,0],[0,1,0],[0,0,1]],dtype=float) U=np.array([[0,0,0],[0,0,0],[0,0,0]],dtype=float) b=np.array([7,-1,0],dtype=float) y=np.array([0,0,0],dtype=float) x=np.array([0,0,0],dtype=float) def LU(A): n=len(A[0]) i=0 while i

相减得01)51()51(e h e h e n n n -==-=-Λ 当 151≤-h 时,4.00≤

常微分方程数值解法

第八章 常微分方程的数值解法 一.内容要点 考虑一阶常微分方程初值问题:?????==0 0)() ,(y x y y x f dx dy 微分方程的数值解:设微分方程的解y (x )的存在区间是[a,b ],在[a,b ]内取一系列节 点a= x 0< x 1<…< x n =b ,其中h k =x k+1-x k ;(一般采用等距节点,h=(b-a)/n 称为步长)。在每个节点x k 求解函数y(x)的近似值:y k ≈y(x k ),这样y 0 , y 1 ,...,y n 称为微分方程的数值解。 用数值方法,求得f(x k )的近似值y k ,再用插值或拟合方法就求得y(x)的近似函数。 (一)常微分方程处置问题解得存在唯一性定理 对于常微分方程初值问题:?????==0 0)() ,(y x y y x f dx dy 如果: (1) 在B y y A x x 00≤-≤≤,的矩形内),(y x f 是一个二元连续函数。 (2) ),(y x f 对于y 满足利普希茨条件,即 2121y y L y x f y x f -≤-),(),(则在C x x 0≤≤上方程?????==0 0)() ,(y x y y x f dx dy 的解存在且唯一,这里C=min((A-x 0),x 0+B/L),L 是利普希茨常数。 定义:任何一个一步方法可以写为),,(h y x h y y k k k 1k Φ+=+,其中),,(h y x k k Φ称为算法的增量函数。 收敛性定理:若一步方法满足: (1)是p 解的. (2) 增量函数),,(h y x k k Φ对于y 满足利普希茨条件. (3) 初始值y 0是精确的。则),()()(p h O x y kh y =-kh =x -x 0,也就是有 0x y y lim k x x kh 0h 0 =--=→)( (一)、主要算法 1.局部截断误差 局部截断误差:当y(x k )是精确解时,由y(x k )按照数值方法计算出来的1~ +k y 的误差y (x k+1)- 1~ +k y 称为局部截断误差。 注意:y k+1和1~ +k y 的区别。因而局部截断误差与误差e k +1=y (x k +1) -y k +1不同。 如果局部截断误差是O (h p+1),我们就说该数值方法具有p 阶精度。

常微分方程的几种数值解法

┊┊┊┊┊┊┊┊┊┊┊┊┊装┊┊┊┊┊订┊┊┊┊┊线┊┊┊┊┊┊┊┊┊┊┊┊┊ 常微分方程的几种数值解法 数学与应用数学肖振华指导教师张秀艳 【摘要】自然界与工程技术中的很多现象,可以归结为微分方程定解问题。其中,常微分方程求解是微分方程的重要基础内容。但是,对于许多的微分方程,往往很难得到甚至不存在精确的解析表达式,这时候,数值解提供了一个很好的解决思路。,针对于此,本文对常微分方程数值解法进行了简单研究,主要讨论了一些常用的数值解法,如欧拉法、改进的欧拉法、Runge—Kutta方法、Adams预估校正法以及勒让德谱方法等,通过具体的算例,结合MA TLAB求解画图,初步给出了一般常微分方程数值解法的求解过程。同时,通过对各种方法的误差分析,让大家对各种方法的特点和适用范围有一个直观的感受。 【关键词】常微分方程数值解法MA TLAB 误差分析 【Abstract】 Many phenomena in nature and engineering can be attributed to the definite solution of the problem for differential equations. Among them, the ordinary differential equation solving is an important foundation for the content of the differential equations. However, many of the differential equations are often difficult to obtain accurate analytical expression .At this time, the numerical solution provides a good idea. For the Numerical Solution of Ordinary Differential Equations in this article, we focuses on some commonly used numerical solution, such as the Euler method, improved Euler method, Runge-Kutta method, Adams predictor corrector method as well as newer spectral methods. Through specific examples, combined with MATLAB solving and drawing, we initially know the solution process of general numerical solution of ordinary differential equations . At the same time, according to the error analysis of various methods , everyone has an intuitive feel of the characteristics and scope of the various methods. 【Keywords】Ordinary Differential Equations Numerical Solution MATLAB error analysis

(整理)常微分方程数值解法

i.常微分方程初值问题数值解法 i.1 常微分方程差分法 考虑常微分方程初值问题:求函数()u t 满足 (,), 0du f t u t T dt =<≤ (i.1a ) 0(0)u u = (i.1b) 其中(,)f t u 是定义在区域G : 0t T ≤≤, u <∞上的函数,0u 和T 是给定的常数。我们假设(,)f t u 对u 满足Lipschitz 条件,即存在常数L 使得 121212(,)(,), [0,]; ,(,)f t u f t u L u u t T u u -≤-?∈∈-∞∞ (i.2) 这一条件保证了(i.1)的解是适定的,即存在,唯一,而且连续依赖于初值0u 。 通常情况下,(i.1)的精确解不可能用简单的解析表达式给出,只能求近似解。本章讨论常微分方程最常用的近似数值解法--差分方法。先来讨论最简单的Euler 法。为此,首先将求解区域[0,]T 离散化为若干个离散点: 0110N N t t t t T -=<<<<=L (i.3) 其中n t hn =,0h >称为步长。 在微积分课程中我们熟知,微商(即导数)是差商的极限。反过来,差商就是微商的近似。在0t t =处,在(i.1a )中用向前差商 10()()u t u t h -代替微商du dt ,便得 10000()()(,())u t u t hf t u t ε=++ 如果忽略误差项0ε,再换个记号,用i u 代替()i u t 便得到 1000(,)u u hf t u -= 一般地,我们有 1Euler (,), 0,1,,1n n n n u u hf t u n N +=+=-L 方法: (i.4) 从(i.1b) 给出的初始值0u 出发,由上式可以依次算出1,,N t t L 上的差分解1,,N u u L 。

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

《计算机数学基础》数值部分第五单元辅导 14 常微分方程的数值解法 一、重点内容 1. 欧拉公式: ),...,,,() ,()(1-210=???+=+=≈01+1+n k kh x x y x hf y y x y k k k k k k 局部截断误差是O (h 2)。 2. 改进欧拉公式: 预报-校正公式: ?? ???++=+=++++)],(),([2),(1111k k k k k k k k k k y x f y x f h y y y x hf y y 校正值 预报值 即 ))],(,(),([2 11k k k k k k k k y x hf y x f y x f h y y +++ =++ 或表成平均的形式: ??? ? ? ???? +21=+=+=1+1+)(),(),(c p k p k k c k k k p y y y y x hf y y y x hf y y 改进欧拉法的局部截断误差是O (h 3) 3. 龙格-库塔法 二阶龙格-库塔法的局部截断误差是O (h 3) 三阶龙格-库塔法的局部截断误差是O (h 4) 四阶龙格?库塔法公式: )22(6 43211κκκκ++++ =+h y y k k 其中 κ1=f (x k ,y k );κ2=f (x n +12h ,y k +21h κ1);κ3=f (x k +12h ,y n +2 1 h κ2);κ4=f (x k +h ,y k +h κ3) 四阶龙格-库塔法的局部截断误差是O (h 5)。 二、实例 例1 用欧拉法解初值问题???1 =060≤≤0--='2)() .(y x xy y y ,取步长h =0.2。计算过程保留 4位小数。 解h =0.2, f (x )=-y -xy 2。首先建立欧拉迭代格式 ) ,,)((.),(210=-420=--=+=2 1+k y x y y hx hy y y x hf y y k k k k k k k k k k k 当k =0,x 1=0.2时,已知x 0=0,y 0=1,有

常微分方程数值解法

常微分方程数值解法 【作用】微分方程建模是数学建模的重要方法,因为许多实际问题的数学描述将导致求解微分方程的定解问题。把形形色色的实际问题化成微分方程的定解问题,大体上可以按以下几步: 1. 根据实际要求确定要研究的量(自变量、未知函数、必要的参数等)并确定坐标系。 2. 找出这些量所满足的基本规律(物理的、几何的、化学的或生物学的等等)。 3. 运用这些规律列出方程和定解条件。基本模型 1. 发射卫星为什么用三级火箭 2. 人口模型 3. 战争模型 4. 放射性废料的处理通常需要求出方程的解来说明实际现象,并加以检验。如果能得到解析形式的解固然是便于分析和应用的,但是我们知道,只有线性常系数微分方程,并且自由项是某些特殊类型的函数时,才可以得到这样的解,而绝大多数变系数方程、非线性方程都是所谓“解不出来” 的于是对于用微分方程解决实际问题来说,数值解法就是一个十分重要的手段。 1. 改进Euler 法: 2. 龙格—库塔( Runge—Kutta )方法: 【源程序】 1. 改进Euler 法: function [x,y]=eulerpro(fun,x0,x1,y0,n);%fun 为函数,(xO, x1)为x 区间,yO 为初始值,n 为子 区间个数 if nargin<5,n=5O;end h=(x1-xO)/n; x(1)=xO;y(1)=yO; for i=1:n x(i+1)=x(i)+h; y1=y(i)+h*feval(fun,x(i),y(i)); y2=y(i)+h*feval(fun,x(i+1),y1); y(i+1)=(y1+y2)/2; end 调用command 窗口 f=i nlin e('-2*y+2*x A2+2*x') [x,y]=eulerpro(f,O,,1,1O) 2 x +2x , (0 < x < , y(0) = 1 求解函数y'=-2y+2 2. 龙格—库塔( Runge—Kutta )方法: [t,y]=solver('F',tspan ,y0) 这里solver为ode45, ode23, ode113,输入参数F是用M文件定义的微分方程y'= f (x, y)右端的函数。tspan=[t0,tfinal]是求解区间,y0是初值。 注:ode45和ode23变步长的,采用Runge-Kutta算法。 ode45表示采用四阶-五阶Runge-Kutta算法,它用4阶方法提供候选解,5阶方法控制误差,是一种自适应步长(变步长)的常微分方程数值解法,其整体截断误差为(△ 口人5解 决的是Nonstiff(非刚性)常微分方程。

相关文档