matlab积分与微分实验报告,常微分matlab实验报告

常微分matlab实验报告

数学与信息科学系实训报告 常微分方程初值问题数值解法及其MATLAB实现 姓名： 学号： 专业：信息与计算科学 年级：2010级 指导教师：朱耀生 完成时间：2013年11月25—2013年12月6 实验目的: (1)研究满足给定方程的可微函数的数值方法 (2)培养matlab编程和上机调试能力 实验基本原理和内容: 根据给定的初始条件，确定常微分方程惟一解的问题叫常微分方程初值问题。大多数实际问题难以求得解析解，必须将微分问题离散化，用数值方法求其近似解。 一阶常微分方程的初值问题的提法是，求出函数，使满足条件 (1) 利用数值方法解问题 (1)时，通常假定解存在且惟一，解函数及右端函数具有所需的光滑程度。数值解法的基本思想是:先取自变量一系列离散点,把微分问题(1)离散化，求出离散问题的数值解,并以此作为微分问题解的近似。例如取步长>0,以剖分区间[,]，令=+,把微分方程离散化成一个差分方程。以()表微分方程初值问题的解,以表差分问题的解，就是近似解的误差，称为全局误差。因此，设计各种离散化模型,求出近似解,估计误差以及研究数值方法的稳定性和收敛性等构成了数值解法的基本内容 1.欧拉法和后退欧拉法 原理: function [x,y]=meuler(df,xspan,y0,h) %用途：改进欧拉公式解常微分方程y =f(x,y), y(x0)=y0 %格式：[xy]=meuler(df,a,b,y0,h) df为函数f(x,y), xspan为求解 %区间[x0,xn], y0为初值y(x0), h为步长, [xy]返回节点和数值解矩阵 x=xspan(1):h:xspan(2); y(1)=y0; for n=1:(length(x)-1) k1=f(df,x(n),y(n)); y(n+1)=y(n)+h*k1; k2=f(df,x(n+1),y(n+1)); y(n+1)=y(n)+h*(k1+k2)/2; end 例题 给定的初值问题y′=-y+x+2，0=h=0.1; y=-1; x=1; for i=1:20 k1=h*oulei_wf(x,y); k2=h*oulei_wf(x+h,y+k1); y=y+0.5*k1+0.5*k2 x=x+h; z=ouleij_q(x) t=y-z end z = 4,z = 3.7000,y = -0.6150 z = 1.4329, z =1.4329, t = -2.0479 z =3.7150, z =3.4435, y = -0.2571 z =1.5012, z =1.5012, t =-1.7583 z =3.4571, z = 3.2114, y =0.0763 z =1.5725, z =1.5725, t =-1.4962 z =3.2237, z = 3.0013, y =0.3876 z =1.6466, z =1.6466, t =-1.2590 z =3.0124,z =2.8112,y =0.6788 z =1.7231,z =1.7231,t =-1.0444 z =2.8212,z = 2.6391,y =0.9518 z =1.8019,z =1.8019,t =-0.8501 z = 2.6482,z =2.4834,y =1.2084 z =1.8827,z =1.8827, t =-0.6743 z =2.4916,z =2.3425, y =1.4501 z =1.9653,z =1.9653, t =-0.5152 实验结果的分析与研究 1.对于欧拉法，步长越小，精度越高，而产生的误差越小。总 对于欧拉法，步长越小，精度越高，而产生的误差越小。 对于欧拉法 体来说，欧拉法的优点是形式简单 计算方便， 是形式简单， 体来说，欧拉法的优点是形式简单，计算方便，缺点是总的运 算精度比较低。 的增大，误差值也越来越大。 算精度比较低。而且随着 x 的增大，误差值也越来越大。根据 欧拉公式的截断误差计算，欧拉法是一阶方法。 欧拉公式的截断误差计算，欧拉法是一阶方法。 2.对于改进欧拉法， 对于改进欧拉法， 其基本特征与欧拉法相似， 也是步长越小， 对于改进欧拉法 其基本特征与欧拉法相似， 也是步长越小， 精度越高，误差越小。优点是精度相对欧拉法来说较高， 精度越高，误差越小。优点是精度相对欧拉法来说较高，截断 误差为 O(h^3)，缺点是比欧拉法计算量大。根据欧拉改进法 ，缺点是比欧拉法计算量大。 公式的截断误差计算，欧拉改进法是二阶方法。 公式的截断误差计算，欧拉改进法是二阶方法。 2．龙格-库塔方法 原理：对于龙哥库塔三阶： 利用Yn+1 = Yn + h/6(K1 + 4K2 + K3 ) K1 = f(Xn, Yn) K2 = f(Xn + 1/2*h, Yn +( h/2)*K1) K3= f(Xn + h, Yn - K1+2*h*k2) 若>> clear all h=0.1; y=-1; x=1; for i=1:25; k1=oulei_wf(x,y); k2=oulei_wf(x+h/2,y+k1*(h/2)); k3=oulei_wf(x+h,y-k1*h+k2*(2*h)); y=y+(k1+4*k2+k3)*(h/6); x=x+h; z=ouleij_q(x); t=abs(y-z); A=[x y z t] end z = 4,z = 3.8500,z = 3.7300, A = 1.1000 -0.6145 1.4329 2.0474 z =3.7145z=3.5788z = 3.4702z =1.5012 A = 1.2000 -0.2562 1.5012 1.7574 z =3.4562z =3.3334z =3.2351z =1.5725 A = 1.3000 0.0776 1.5725 1.4950 z = 3.2224z =3.1113z =3.0224z =1.6466 A = 1.4000 0.3891 1.6466 1.2575 z = 3.0109z = 2.9104z =2.8299z =1.7231 A =1.5000 0.6804 1.7231 1.0427 z =2.8196z =2.7286z =2.6558z =1.8019 A = 1.6000 0.9536 1.8019 0.8483