微分方程数值解试验报告.doc

中国矿业大学（北京）理学院微分方程数值解实验报告实验名称： 常微分方程数值解法 学号： XXX 专业年级： 信息与计算科学11级 姓名： XXX 实验时间： 2013年3月 成绩 ： 一、实验目的，内容 二、相关背景知识介绍 三、代码四、数值结果与分析 六、计算结果的分析 1、 实验目的，内容 掌握一阶常微分方程的初值问题常用数值解法 编写程序，利用Euler法、改进Euler法、Adams外插内插法及Runge-Kutta法解决相关问题。2、 相关背景知识介绍 研究一阶常微分方程的初值问题的数值解 1、Euler法和改进Euler法 （1）Euler法 用差商近似导数 问题转化为（2）改进Euler法 （3）线性多步法（Adams外插内插法） （4）Runge-Kutta法 用f在一些点的值表示（tn,un,h）,使单步局部截断误差的阶和Taylor展开法相等。将初值问题写成积分形式：3、 代码（Matlab） 题目：用Euler法、改进Euler法、Adams外插内插法及Runge-Kutta法解决下面一阶常微分方程初值问题。代码:方程： function y=f(t,u) y=2*t*u/(1+t2); endEuler法: function y=euler(n,u0,x0) t=1; u=zeros(n+1,1); h=t/n; u(1)=u0; x(1)=x0; for i=1:n x(i+1)=x0+h*i; u(i+1)=u(i)+f(x(i),u(i)*h; end Plot(x,u,r*) end改进Euler法: function y=imeuler(n,u0,x0) err=10-6; t=1; u=zeros(n+1,1); x=zeros(n+1,1); h=t/n; u(1)=u0; x(1)=x0; for i=1:n x(i+1)=x0+h*i; u(i+1)=u(i); ut=30; while abs(u(i+1)-ut)err ut=u(i+1); u(i+1)=u(i)+(f(x(i),u(i)+f(x(i+1),u(i+1)*0.5*h; end end plot(x,u) endAdams内插法：k=3function y=iadams(n,u0,x0) t=1; u=zeros(n+1,1); h=t/n; u(1)=u0; u(2)=u0; u(3)=u0; x(1)=x0; x(2)=x0+h; x(3)=x0+2*h; for i=3:n x(i+1)=x0+h*i; u(i+1)=u(i)+h*(9*f(x(i+1),u(i+1)+19*f(x(i),u(i)-5*f(x(i-1),u(i-1)+f(x(i-2),u(i-2)/24; end plot(x,u) endRunge-Kutta法：（本文采用heun法）function y=heun(n,u0,x0) t=1; u=zeros(n+1,1); h=t/n; u(1)=u0; x(1)=x0; for i=1:n x(i+1)=x0+h*i; k1=f(x(i),u(i); k2=f(x(i)+0.3333*h,u(i)+0.33333*h*k1); k3=f(x(i)+0.666663*h,u(i)+0.666666*h*k2); u(i+1)=u(i)+0.25*h*(k1+3*k3); end plot(x,u) end 四、数值结果与分析 对于所给方程，各个方法均给出了良好的解，其实heun方法最好10步就出了 的结果。而euler法10步 Euler法大致形象出来了，但是后边误差太大。Adams内插法一开始浪费了三个点的 10