实验报告七常微分方程初值问题数值解法
浙江高校城市学院试验报告 课程名称 数值计算方法 试验项目名称 常微分方程初值问题的数值解法 试验成果 指导老师(签名 ) 日期 2015/12/16 一. 试验目的和要求 1. 用Matlab软件驾驭求微分方程数值解的欧拉方法和龙格-库塔方法; 2. 通过实例学习用微分方程模型解决简化的实际问题。 二. 试验内容和原理 编程题2-1要求写出Matlab源程序(m文件),并有适当的注释语句;分析应用题2-2,2-3,2-4,2-5要求将问题的分析过程、Matlab源程序和运行结果和结果的说明、算法的分析写在试验报告上。 2-1 编程 编写用向前欧拉公式和改进欧拉公式求微分方程数值解的Matlab程序,问题如下: 在区间内个等距点处,靠近下列初值问题的解,并对程序的每一句添上注释语句。 Euler法 y=euler(a,b,n,y0,f,f1,b1) 改进Euler法 y=eulerpro(a,b,n,y0,f,f1,b1) 2-2 分析应用题 假设等分区间数,用欧拉法和改进欧拉法在区间内求解初值问题 并作出解的曲线图形,同时将方程的解析解也画在同一张图上,并作比较,分析这两种方法的精度。 2-3 分析应用题 用以下三种不同的方法求下述微分方程的数值解,取 画出解的图形,与精确值比较并进行分析。 1)欧拉法; 2)改进欧拉法; 3)龙格-库塔方法; 2-4 分析应用题 考虑一个涉及到社会上别出心裁的人的繁衍问题模型。假设在时刻(单位为年),社会上有人口人,又假设全部别出心裁的人与别的别出心裁的人结婚后所生后代也是别出心裁的人。而固定比例为的全部其他的后代也是别出心裁的人。假如对全部人来说诞生率假定为常数,又假如一般的人和别出心裁的人的婚配是随意的,则此问题可以用微分方程表示为: 其中变量表示在时刻社会上别出心裁的人的比例,表示在时刻人口中别出心裁的人的数量。 1)假定和,当步长为年时,求从到解的近似值,并作出近似解的曲线图形。 2)精确求出微分方程的解,并将你当时在分题(b)中得到的结果与此时的精确值进行比较。 【MATLAB相关函数】 n 求微分方程的解析解及其数值的代入 dsolve(‘egn1’, ‘egn2’, ‘’) subs (expr, {x,y,…}, {x1,y1,…} ) 其中‘egn’表示第个方程,‘’表示微分方程中的自变量,默认时自变量为。 subs吩咐中的expr、x、y为符合型表达式,x、y分别用数值x1、x2代入。 >> syms x y z >> subs( x+y+z ,{x,y,z},{1,2,3}) ans = 6 >> syms x >> subs( x^2 ,x,2) ans = 4 s=dsolve(‘’, ‘’, ‘’) ans = >> syms x >> subs(s,x,2) ans = -0.3721 n 右端函数的自动生成 f= inline(‘expr’, ’var1’, ‘var2’,……) 其中’expr’表示函数的表达式,’var1’, ‘var2’ 表示函数表达式中的变量,运行该函数,生成一个新的函数表达式为f (var1, var2, ……)。 >> f=inline( x+3*y , x , y ) f = Inline function: f(x,y) = x+3*y >> f(2,3) ans = 11 n 4,5阶龙格-库塔方法求解微分方程数值解 [t,x]=ode45(f,ts,x0,options) 其中f是由待解方程写成的m文件名;x0为函数的初值;t,x分别为输出的自变量和函数值(列向量),t的步长是程序依据误差限自动选定的。若ts=[t0,t1,t2,…,tf],则输出在自变量指定值,等步长时用ts=t0:k:tf,输出在等分点;options用于设定误差限(可以缺省,缺省时设定为相对误差,肯定误差),程序为: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,off off 例:解微分方程 在吩咐窗口执行 = (‘’, ‘’, ‘’); ; ans = 0 1.0000 0.0502 1.0490 0.1005 1.0959 0.1507 1.1408 3.8507 2.9503 3.9005 2.9672 3.9502 2.9839 4.0000 3.0006 plot(,,‘o-’,) %解函数图形表示 %不用输出变量,则干脆输出图形 ; ans = 0 1.0000 1.0000 1.7321 2.0000 2.2361 3.0000 2.6458 4.0000 3.0006 三. 操作方法与试验步骤(包括试验数据记录和处理) 2-1 编程 编写用向前欧拉公式和改进欧拉公式求微分方程数值解的Matlab程序,问题如下: 在区间内个等距点处,靠近下列初值问题的解,并对程序的每一句添上注释语句。 Euler法 y=euler(