微分方程数值解法实验报告
常微分方程数值解实验报告
常微分方程数值解实验报告学院：数学与信息科学专业：信息与计算科学姓名：郑思义学号：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进行计算，并与真解作比较。
实验报告——常微分方程的数值解法
实验报告实验项目名称常微分方程的数值解法实验室数学实验室所属课程名称微分方程数值解实验类型上机实验实验日期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 法的精度更高。
常微分方程数值解实验报告
常微分方程数值解实验报告实验报告：常微分方程数值解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尝试不同的步长，观察相应的数值解的变化。
大学数学实验四_常微分方程数值解
2 小型火箭初始质量为 1400kg，其中包括 1080kg 燃料。火箭竖直向上发射时，燃料燃烧率 为 18kg/s，由此产生 32000N 推力，火箭引擎在燃料用尽时关闭。设火箭上升时空气阻力正 比于速度的平方，比例系数为 0.4kg/m，求引擎关闭瞬间火箭的高度、速度、加速度，及火 箭到达最高点时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。
6
输出结果为
故，当 t =60s 时，火箭的速度、高度、加速度为
v (m/s)
向前欧拉公式
266.94
经典龙格-库塔公式
266.94
ode45(@rocket,t,v)
266.93
h (m) 12131 12131 12131
从上表可以看出，用三种不同方法求得的结果相差不大。
60s 后，火箭作竖直上抛运动，加速度为
4
输出的结果为
输出的火箭的 h – t 图线如下：
经典龙格-库塔公式： n=100000; h=60/n; t=0:h:60; v(1)=0; z(1)=0; g(1)=90/7; for i=1:n
a = v(i); s = z(i); x1 = h*i; k1 = (32000-0.4*a.^2)/(1400-18*x1)-10; x1 = x1 + h/2; b = a + h/2*k1; k2 = (32000-0.4*b.^2)/(1400-18*x1)-10; c = a + h/2*k2; k3 = (32000-0.4*c.^2)/(1400-18*x1)-10; x1 = x1 + h/2;
n = 500 267.0539 n = 100000 266.9432
微分方程数值解法实验报告
微分方程数值解法实验报告班级：姓名：学号：日期：一、实验目的1、熟悉微分方程（组）数值解的Euler算法，改进的Euler算法和Runge-Kutta算法，利用matlab软件实现微分方程数值解法来求解具体试题；2、虽然求解常微分方程有各种各样的解析解，但解析方法只能用来求解一些特殊类型的方程，通常它们无法求出解析解，而需要数值方法来近似求解。
因此产生了常微分方程初值问题的数值计算方法，常微分方程数值解法是通过计算机便捷的求解近似值。
二、基本理论及背景1、在数值求解常微分方程中，主要有有限差分计算和有限元计算两大类方法，其中在有限差分计算方法中有一类方法称为龙格－库塔(Runge_Kutta)方法。
四阶的龙格－库塔方法为最佳的计算格式。
2、参考三中的代码，分别用Euler算法，改进的Euler算法和Runge-Kutta 算法实现微分方程（组）的数值求解，完成下列题目：三、算法设计及实现1、算法设计，通过Euler算法，改进的Euler算法和Runge-Kutta三种算法来实现微分方程（组）的数值求解；2、程序文件及功能清单：(1) Euler Method:function [x,y]=EulerDSolve(f,ab,y0,h)x=(ab(1):h:ab(2))';n=length(y0);y=zeros(length(x),n);y(1,:)=y0';for k=2:length(x)y(k,:)=y(k-1,:)+h*feval(f,x(k-1),y(k-1,:)')';end；(2) Improved Euler Method:function [x,y]=MEulerDSolve(f,ab,y0,h)x=(ab(1):h:ab(2))';n=length(x);y=zeros(n,length(y0));y(1,:)=y0';for k=2:nyp=y(k-1,:)+h*feval(f,x(k-1),y(k-1,:));yc=y(k-1,:)+h*feval(f,x(k),yp);y(k,:)=(yp+yc)/2;end(3) Runge-Kutta Method:function [TOut,YOut]=Runge_Kutta(f,ab,y0,h)TOut=(ab(1):h:ab(2))';n=length(TOut);YOut=zeros(n,length(y0));YOut(1,:)=y0';for k=2:nx=TOut(k-1); y=YOut(k-1,:)';K1=feval(f,x,y);K2=feval(f,x+h/2,y+K1*h/2);K3=feval(f,x+h/2,y+K2*h/2);K4=feval(f,x+h,y+K3*h);YOut(k,:)=(y+(K1+2*K2+2*K3+K4)*h/6)';end四、实验步骤1、打开MATLAB软件，新建 *.m文件，在m文件的窗口中编辑Euler算法的函数程序，另建一m文件，编辑自己改进的Euler算法的函数程序，再新建一m文件，在窗口中编辑Runge-Kutta算法的函数程序，并全部保存在指定的文件夹下；2、将MATLAB软件的工作页面的工具栏下的目标文件指向指定的文件夹；3、分别调用上述三种算法的函数，实现微分方程（组）的数值求解完成给定的实验题目；4、输出结果和初步分析说明（见附页）。
微分方程数值解法实验报告
微分方程数值解法实验报告2篇微分方程数值解法实验报告（一）在实际科学与工程问题中，我们经常会遇到微分方程的求解。
然而，许多微分方程往往没有解析解，这就需要我们利用数值方法来获得近似解。
本篇实验报告将介绍两种常见的微分方程数值解法：欧拉方法和改进的欧拉方法。
一、欧拉方法欧拉方法是最简单的微分方程数值解法之一。
其基本原理为离散化微分方程，将微分方程中的导数用差商代替。
设要求解的微分方程为dy/dx = f(x, y)，步长为h，则可用以下公式进行递推计算：y_{n+1} = y_n +hf(x_n, y_n)二、算法实现为了对欧拉方法进行数值实验，我们以一阶线性常微分方程为例：dy/dx = x - y, y(0) = 1步骤如下：（1）选择合适的步长h和求解区间[a, b]，这里我们取h=0.1，[a, b] = [0, 1]；（2）初始化y_0 = 1；（3）利用欧拉方法递推计算y_{n+1} = y_n + 0.1(x_n - y_n)；（4）重复步骤（3），直到x_n = 1。
三、实验结果与讨论我们通过上述步骤得到了在[0, 1]上的近似解y(x)。
下图展示了欧拉方法求解的结果。
从图中可以看出，欧拉方法得到的近似解与精确解有一定的偏差。
这是因为欧拉方法只是通过递推计算得到的近似解，并没有考虑到更高阶的误差。
所以在需要高精度解时，欧拉方法并不理想。
四、改进的欧拉方法针对欧拉方法的不足，我们可以考虑使用改进的欧拉方法（也称为改进的欧拉-柯西方法）。
它是通过利用前后两个步长欧拉方法得到的结果作为差商的中间项，从而提高了求解精度。
一阶线性常微分方程的改进欧拉方法可以表示为：y_{n+1} = y_n + h(f(x_n, y_n) + f(x_{n+1}, y_n + hf(x_n,y_n)))/2五、算法实现与结果展示将改进的欧拉方法应用于前述的一阶线性常微分方程，我们同样选择h=0.1，[a, b] = [0, 1]。
微分方程数值解实验报告
微分方程数值解实验报告实验目的：掌握微分方程数值解的基本方法，能够利用计算机编程求解微分方程。
实验原理：微分方程是自然科学与工程技术中常见的数学模型，它描述了变量之间的关系及其随时间、空间的变化规律。
解微分方程是研究和应用微分方程的基础，但有很多微分方程无法找到解析解，只能通过数值方法进行求解。
本实验采用欧拉方法和改进的欧拉方法求解微分方程的初值问题：$$\begin{cases}\frac{dy}{dt}=f(t,y)\\y(t_0)=y_0\end{cases}$$其中，$f(t,y)$是给定的函数，$y(t_0)=y_0$是已知的初值条件。
欧拉方法是最基本的数值解法，其步骤如下：1.给定$t_0$和$y_0$2.计算$t_{i+1}=t_i+h$，其中$h$是步长3. 计算$y_{i+1}=y_i+hf(t_i,y_i)$4.重复步骤2、3直到达到终止条件改进的欧拉方法是对欧拉方法进行改进，通过利用函数$y(t)$在$t+\frac{1}{2}h$处的斜率来更准确地估计$y_{i+1}$，其步骤如下：1.给定$t_0$和$y_0$2.计算$t_{i+1}=t_i+h$，其中$h$是步长3. 计算$y_*=y_i+\frac{1}{2}hf(t_i,y_i)$4. 计算$y_{i+1}=y_i+hf(t_i+\frac{1}{2}h,y_*)$5.重复步骤2、3、4直到达到终止条件实验步骤：1.编写程序实现欧拉方法和改进的欧拉方法2.给定微分方程和初值条件3.设置步长和终止条件4.利用欧拉方法和改进的欧拉方法求解微分方程5.比较不同步长下的数值解与解析解的误差6.绘制误差-步长曲线，分析数值解的精度和收敛性实验结果：以一阶常微分方程$y'=3ty+t$为例，给定初值$y(0)=1$，取步长$h=0.1$进行数值求解。
利用欧拉方法求解微分方程得到的数值解如下：\begin{array}{cccc}t & y_{\text{exact}} & y_{\text{Euler}} & \text{误差} \\ \hline0.0&1.000&1.000&0.000\\0.1&1.035&1.030&0.005\\0.2&1.104&1.108&0.004\\0.3&1.212&1.217&0.005\\0.4&1.360&1.364&0.004\\0.5&1.554&1.559&0.005\\0.6&1.805&1.810&0.005\\0.7&2.131&2.136&0.005\\0.8&2.554&2.560&0.006\\0.9&3.102&3.107&0.006\\1.0&3.807&3.812&0.005\\\end{array}利用改进的欧拉方法求解微分方程得到的数值解如下：\begin{array}{cccc}t & y_{\text{exact}} & y_{\text{Improved Euler}} & \text{误差} \\\hline0.0&1.000&1.000&0.000\\0.1&1.035&1.035&0.000\\0.2&1.104&1.103&0.001\\0.3&1.212&1.211&0.001\\0.4&1.360&1.358&0.002\\0.5&1.554&1.552&0.002\\0.6&1.805&1.802&0.003\\0.7&2.131&2.126&0.005\\0.8&2.554&2.545&0.009\\0.9&3.102&3.086&0.015\\1.0&3.807&3.774&0.032\\\end{array}误差-步长曲线如下：实验结论：通过对比欧拉方法和改进的欧拉方法的数值解与解析解的误差，可以发现改进的欧拉方法具有更高的精度和收敛性。
微分方程数值方法实验报告
微分方程数值方法实验报告微分方程数值方法实验报告一、实验目的1、了解Euler法及梯形法的基本原理，能用其解决常微分方程初值问题，并把计算结果用图形表示出来。
2、理解4阶RK法基本计算步骤，能编程实现算法并解决相关常微分方程初值问题。
3、了解MATLAB主要功能和基本特征，熟悉MATLAB操作环境。
掌握MATLAB常用函数的使用以及图形处理。
二、实验题目对于初值问题u’=u,u(0)=1,在区间[0,1]上，用Euler法，梯形法及RK方法进行计算，取步长h=0.1,0.2,0.5，试比较(1)用同样步长，三种方法中哪一个精度最好;(2)对同一种方法一不同步长进行计算，哪一个结果最好。
三、实验内容1、步长为h=0.1时，用三种方法计算题目1)、MATLAB程序Euler法:>> a=0;b=1;h=0.1;N=(b-a)/h;y=zeros(N+1,1);y(1)=1;x=a:h:b;>> for i=2:N+1y(i)=y(i-1)+h*y(i-1);end求得:y = (1 1.1 1.21 1.331 1.4641 1.6105 1.7716 1.9487 2.1436 2.3579 2.5937)’梯形法:>> z=zeros(N+1,1);>> z(1)=1;>> for i=2:N+1v1=z(i-1);t=z(i-1)+h*v1;v2=t;z(i)=z(i-1)+h/2*(v1+v2);end1求得:z = (1 1.105 1.2205 1.3476 1.4873 1.641 1.8101 1.99622.2008 2.4258 2.6734)’RK法:>> w=zeros(N+1,1);>> w(1)=1;>> for i=2:N+1K1=w(i-1);K2=w(i-1)+K1*h/2;K3=w(i-1)+K2*h/2;K4=w(i-1)+K3*h;w(i)=w(i-1)+h*(K1+2*K2+2*K3+K4)/6;end求得:w =(1 1.1052 1.2214 1.3499 1.4918 1.6487 1.8221 2.01382.2255 2.4596 2.7183)’>> plot(x,y)>> hold on>> plot(x,z,':')>> plot(x,w,'--')>> plot(x,exp(x),'*')>> title('相同步长下三种方法与准确解的对比')>> legend('Euler法','梯形法','四阶RK法','准确解') 2)、图形对比相同步长下三种方法与准确解的对比2.82.62.42.221.81.6Euler法1.4梯形法四阶RK法1.2准确解1 00.10.20.30.40.50.60.70.80.9123)、结果分析由图得出三种方法中四阶RK法经确度最高，梯形法次之，Euler法精确度最差。
实验报告七常微分方程初值问题的数值解法
浙江大学城市学院实验报告课程名称数值计算方法实验项目名称常微分方程初值问题的数值解法 实验成绩指导老师签名日期2015/12/16 一.实验目的和要求1． 用Matlab 软件掌握求微分方程数值解的欧拉方法和龙格－库塔方法； 2． 通过实例学习用微分方程模型解决简化的实际问题;二.实验内容和原理编程题2-1要求写出Matlab 源程序m 文件,并有适当的注释语句；分析应用题2-2,2-3,2-4,2-5要求将问题的分析过程、Matlab 源程序和运行结果和结果的解释、算法的分析写在实验报告上; 2-1 编程编写用向前欧拉公式和改进欧拉公式求微分方程数值解的Matlab 程序,问题如下：在区间[],a b 内(1)N +个等距点处,逼近下列初值问题的解,并对程序的每一句添上注释语句; Euler 法y=eulera,b,n,y0,f,f1,b1改进Euler 法y=eulerproa,b,n,y0,f,f1,b1 2-2 分析应用题假设等分区间数100n =,用欧拉法和改进欧拉法在区间[0,10]t ∈内求解初值问题()()20(0)10y t y t y '=-⎧⎨=⎩并作出解的曲线图形,同时将方程的解析解也画在同一张图上,并作比较,分析这两种方法的精度; 2-3 分析应用题用以下三种不同的方法求下述微分方程的数值解,取10h = 画出解的图形,与精确值比较并进行分析; 1欧拉法； 2改进欧拉法； 3龙格－库塔方法；2-4 分析应用题考虑一个涉及到社会上与众不同的人的繁衍问题模型;假设在时刻t 单位为年,社会上有人口()x t 人,又假设所有与众不同的人与别的与众不同的人结婚后所生后代也是与众不同的人;而固定比例为r 的所有其他的后代也是与众不同的人;如果对所有人来说出生率假定为常数b ,又如果普通的人和与众不同的人的婚配是任意的,则此问题可以用微分方程表示为：其中变量()()()i p t x t x t =表示在时刻t 社会上与众不同的人的比例,()i x t 表示在时刻t 人口中与众不同的人的数量;1假定(0)0.01,0.02p b ==和0.1r =,当步长为1h =年时,求从0t =到50t =解()p t 的近似值,并作出近似解的曲线图形;2精确求出微分方程的解()p t ,并将你当50t =时在分题b 中得到的结果与此时的精确值进行比较; MATLAB 相关函数求微分方程的解析解及其数值的代入dsolve‘egn1’,‘egn2’,‘x ’ subsexpr,{x,y,…},{x1,y1,…}其中‘egn i ’表示第i 个方程,‘x ’表示微分方程中的自变量,默认时自变量为t ; subs 命令中的expr 、x 、y 为符合型表达式,x 、y 分别用数值x1、x2代入; >>symsxyz>>subs'x+y+z',{x,y,z},{1,2,3} ans= 6>>symsx>>subs'x^2',x,2 ans= 4>>s=dsolve‘12Dy y ∧=+’,‘(0)1y =’,‘x ’ ans= >>symsx >>subss,x,2 ans=右端函数(,)f x y 的自动生成f=inline ‘expr ’,’var1’,‘var2’,……其中’expr ’表示函数的表达式,’var1’,‘var2’表示函数表达式中的变量,运行该函数,生成一个新的函数表达式为fvar1,var2,……; >>f=inline'x+3y','x','y' f=Inlinefunction: fx,y=x+3y >>f2,3 ans= 114,5阶龙格－库塔方法求解微分方程数值解t,x=ode45f,ts,x0,options其中f 是由待解方程写成的m 文件名；x0为函数的初值；t,x 分别为输出的自变量和函数值列向量,t的步长是程序根据误差限自动选定的;若ts=t0,t1,t2,…,tf,则输出在自变量指定值,等步长时用ts=t0:k:tf,输出在等分点；options 用于设定误差限可以缺省,缺省时设定为相对误差310-,绝对误差610-,程序为：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,offoff例：在命令窗口执行>>odefun =inline ‘2*y t y -’,‘t ’,‘y ’；>>[],45(,[0,4],1)t y ode odefun =；ans=>>t y ‘o-’,%解函数图形表示>>45(,[0,4],1)ode odefun %不用输出变量,则直接输出图形 >>[],45(,0:4,1)t y ode odefun =；[],t yans=三.操作方法与实验步骤包括实验数据记录和处理2-1编程编写用向前欧拉公式和改进欧拉公式求微分方程数值解的Matlab 程序,问题如下：在区间[],a b 内(1)N +个等距点处,逼近下列初值问题的解,并对程序的每一句添上注释语句; Euler 法y=eulera,b,n,y0,f,f1,b1改进Euler 法y=eulerproa,b,n,y0,f,f1,b1Euler 法y=eulera,b,n,y0,f,f1,b1 y=zeros1,n+1; y1=y0; h=b-a/n; x=a:h:b; fori=1:n; yi+1=yi+hfxi,yi; end plotx,y holdon%求微分方程的精确解 x1=linspacea,b,100; '精确解为' s=dsolvef1,b1,'x' symsxy1=zeros1,100; for i=1:100y1i=subss,x,x1i; endplotx1,y1,'r'title'红色代表精确解'改进Euler 法y=eulerproa,b,n,y0,f,f1,b1 %求微分方程的数值解 y=zeros1,n+1; y1=y0; h=b-a/n; x=a:h:b; fori=1:n; T1=fxi,yi; T2=fxi+1,yi+hT1; yi+1=yi+h/2T1+T2; end plotx,y holdon%求微分方程的精确解 x1=linspacea,b,100; '精确解为' s=dsolvef1,b1,'x' symsxy1=zeros1,100; fori=1:100 y1i=subss,x,x1i; endplotx1,y1,'r'title'红色代表精确解' 2-2分析应用题假设等分区间数100n =,用欧拉法和改进欧拉法在区间[0,10]t ∈内求解初值问题()()20(0)10y t y t y '=-⎧⎨=⎩并作出解的曲线图形,同时将方程的解析解也画在同一张图上,并作比较,分析这两种方法的精度;1向前欧拉法>>euler0,10,100,10,inline'y-20','x','y','Dy=y-20','y0=10' ans= 精确解为 s= 20-10expx ans= +005Columns1through8（2）改进欧拉法>>eulerpro0,10,100,10,inline'y-20','x','y','Dy=y-20','y0=10' ans= 精确解为 s= 20-10expx ans= +005Columns1through8改进欧拉法的精度比向前欧拉法更高; 2-3分析应用题用以下三种不同的方法求下述微分方程的数值解,取10h = 画出解的图形,与精确值比较并进行分析; 1欧拉法； 2改进欧拉法；2-4分析应用题考虑一个涉及到社会上与众不同的人的繁衍问题模型;假设在时刻t 单位为年,社会上有人口()x t 人,又假设所有与众不同的人与别的与众不同的人结婚后所生后代也是与众不同的人;而固定比例为r 的所有其他的后代也是与众不同的人;如果对所有人来说出生率假定为常数b ,又如果普通的人和与众不同的人的婚配是任意的,则此问题可以用微分方程表示为：其中变量()()()i p t x t x t =表示在时刻t 社会上与众不同的人的比例,()i x t 表示在时刻t 人口中与众不同的人的数量;1假定(0)0.01,0.02p b ==和0.1r =,当步长为1h =年时,求从0t =到50t =解()p t 的近似值,并作出近似解的曲线图形;2精确求出微分方程的解()p t ,并将你当50t =时在分题b 中得到的结果与此时的精确值进行比较;1>>euler0,50,50,,inline'','t','p','Dp=','p0= 1' ans= 精确解为 s=1-99/100expx/500 ans=Columns1through82>>dsolve'Dp=','p0=','t' ans=1-99/100expt/500 >>1-99/100exp ans=与欧拉法求得的精确值差0,0001四.实验结果与分析。
微分方程数值解实验报告
微分方程数值解实验报告微分方程数值解实验报告一、引言微分方程是数学中一类重要的方程，广泛应用于各个科学领域。
在实际问题中，往往难以得到微分方程的解析解，因此需要借助数值方法来求解。
本实验旨在通过数值解法，探索微分方程的数值解及其应用。
二、数值解法介绍常用的微分方程数值解法有欧拉法、改进欧拉法、四阶龙格-库塔法等。
在本实验中，我们将采用改进欧拉法进行数值解的求取。
改进欧拉法是一种一阶的显式迭代法，其基本思想是将微分方程的导数用差商来近似表示，并通过迭代逼近真实解。
具体迭代公式如下：\[y_{n+1} = y_n + h \cdot f(x_n + \frac{h}{2}, y_n + \frac{h}{2} \cdot f(x_n, y_n))\]其中，\(y_n\)表示第n步的近似解，\(h\)表示步长，\(f(x_n, y_n)\)表示微分方程的导数。
三、实验步骤1. 确定微分方程及初始条件在本实验中，我们选择经典的一阶常微分方程：\[y' = -2xy\]并给定初始条件\(y(0) = 1\)。
2. 设置步长和迭代次数为了得到较为准确的数值解，我们需要合理选择步长和迭代次数。
在本实验中，我们将步长设置为0.1，迭代次数为10。
3. 迭代计算数值解根据改进欧拉法的迭代公式，我们可以通过编写计算程序来求解微分方程的数值解。
具体计算过程如下：- 初始化：设定初始条件\(y_0 = 1\)，并给定步长\(h = 0.1\)。
- 迭代计算：使用改进欧拉法的迭代公式，通过循环计算得到\(y_1, y_2, ...,y_{10}\)。
- 输出结果：将计算得到的数值解输出，并进行可视化展示。
四、实验结果与分析通过以上步骤，我们得到了微分方程的数值解\(y_1, y_2, ..., y_{10}\)。
将这些数值解进行可视化展示，可以更直观地观察到解的变化趋势。
根据实验结果，我们可以发现随着迭代次数的增加，数值解逐渐逼近了真实解。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
微分方程数值解法实验报告姓名： 班级： 学号：一：问题描述求解边值问题：()2(sin cos cos sin (0,1)(0,1)0,(,)x y u e x y x y G u x y G ππππππ+⎧⎫∆=+⎪⎪∈=⨯⎨⎬⎪⎪=∈∂⎩⎭（x,y) 其精确解为)sin()sin(),()(y x e y x u y x πππ+=问题一：取步长h=k=1/64,1/128,作五点差分格式，用Jacobi 迭代法，Gauss_Seidel 迭代法，SOR 迭代法（w=1.45）。
求解差分方程，以前后两次重合到小数点后四位的迭代值作为解的近似值，比较三种解法的迭代次数以及差分解)128/1,64/1)(,(=h y x u h 与精确解的精度。
问题二：取步长h=k=1/64,1/128,作五点差分格式，用单参数和双参数PR 法解差分方程，近似到小数点后四位。
与SOR 法比较精度和迭代步数。
问题三：取步长h=k=1/64,1/128,作五点差分格式，用共轭梯度法和预处理共轭梯度法解差分方程，近似到小数点后四位。
与SOR法与PR 法比较精度和迭代步数。
二．实验目的：分别使用五点差分法（Jacobi 迭代，Gauss_Seidel 迭代，SOR迭代），PR 交替隐式差分法（单参数，双参数），共轭梯度法，预共轭梯度法分别求椭圆方程的数值解。
三．实验原理：(1) Jacobi 迭代法设线性方程组(1)的系数矩阵A 可逆且主对角元素均不为零,令 并将A 分解成(2) 从而(1)可写成令其中. (3) 以为迭代矩阵的迭代法(公式)(4)称为雅可比(Jacobi)迭代法(公式),用向量的分量来表示,(4)为(5) 其中为初始向量. (2) Guass-Seidel 迭代法由雅可比迭代公式可知,在迭代的每一步计算过程中是用的全部分量来计算的所有分量,显然在计算第i 个分量时,已经b Ax =nn a ,...,a ,a 2211()nn a ,...,a ,a diag D 2211=()D D A A +-=()b x A D Dx +-=11f x B x +=b D f ,A D I B 1111--=-=1B ()()111f x B x k k +=+⎩⎨⎧[],...,,k ,n ,...,i x a b a x n ij j )k (j j i i ii )k (i 21021111==∑-=≠=+()()()()()Tn x ,...x ,x x 002010=()k x ()1+k x ()1+k i x计算出的最新分量没有被利用，从直观上看,最新计算出的分量可能比旧的分量要好些.因此，对这些最新计算出来的第次近似的分量加以利用,就得到所谓解方程组的高斯—塞德（Gauss-Seidel ）迭代法.把矩阵A 分解成(6)其中,分别为的主对角元除外的下三角和上三角部分,于是,方程组(1)便可以写成 即其中(7)以为迭代矩阵构成的迭代法(公式)(8)称为高斯—塞德尔迭代法(公式),用 量表示的形式为(3) SOR 迭代(4) 交替方向迭代法（PR 法）迭代格式为：()()1111+-+k i k x ,...,x 1+k ()1+k x ()1+k j x U L D A --=()nn a ,...,a ,a diag D 2211=U ,L --A ()b Ux x L D +=-22f x B x +=()()b L D f ,U L D B 1212---=-=2B ()()221f x B x k k +=+⎩⎨⎧[],...,,k ,n ,,i x a x a b a x i j n i j )k (j ij )k (j ij i ii )k (i 21021111111==∑∑--=-=+=++Λ))1(()(1D R L D T ωωω-+-=-b )(1--=L D d ωωhu πμωcos )11/(22opt =-+=2121,,1,1,1,,122L L L L u u u L u u u j i j i j i j i j i j i +==+-=+-+-+-对于单参数PR 法，对于多参数，(5) 共轭梯度法 算法步骤如下： ［预置步］任意，计算，并令取：指定算法终止常数，置，进入主步；［主步］ （１）如果，终止算法，输出；否则下行；（２）计算：（３）计算：（４）置，转入（１）．(6) 预共轭梯度法b uL I uL I b u L I uL I k k k k k k k k k k ττττττ+-=++-=++++211122211)()()()(hh optπτsin 22=2sin a ....2,1)11(421k 221h k a h k πρρτ==+-=--其中［预置步］任意，计算，并令取：指定算法终止常数，置，进入主步；［主步］（１）计算：，（２）如果，转入（3）．否则，终止算法，输出计算结果（３）计算：（４）置，转入（1）注：在算法［主步］中，引入变量，及，可以简化计算。
四．程序设计(MATLAB实现)Jacobi迭代法function[u_1,m_1]=Jacobi_Solve(A,b,n,err)D=diag(A);D=diag(D);L=-tril(A,-1);R=-triu(A,1);B=D\(L+R);g=(D\b);u_1_0=zeros(n-1,n-1);%初始迭代值u_1_0=u_1_0(:);flag=1;while flagu_1=B*u_1_0+g;if norm(u_1-u_1_0,inf)<errflag=0;endu_1_0=u_1;m_1=m_1+1;%迭代次数值enduu=zeros(n+1);for mm=1:n-1for nn=1:n-1uu(nn+1,mm+1)=u_1((mm-1)*(n-1)+nn,1);endend%Jacobi迭代差分解图像x=[0:1/n:1];y=[0:1/n:1];mesh(x,y,uu);title('Jacobi迭代差分解图像')Gauss_Seidel迭代法function[u_2,m_2]=Gauss_Seidel_Solve(A,b,n,err) if nargin==2err=1e-6;endD=diag(A);D=diag(D);L=-tril(A,-1);U=-triu(A,1);R=(D-L)\U;g=(D-L)\b;m_2=0;u_2_0=zeros(n-1,n-1);%初始迭代值u_2_0=u_2_0(:);flag=1;u_2=R*u_2_0+g;if norm(u_2-u_2_0,inf)<errflag=0;endu_2_0=u_2;m_2=m_2+1;%迭代次数值enduu=zeros(n+1);for mm=1:n-1for nn=1:n-1uu(nn+1,mm+1)=u_2((mm-1)*(n-1)+nn,1);endend%Gauss-Seidel迭代差分解图像x=[0:1/n:1];y=[0:1/n:1];figure(3)mesh(x,y,uu); title('Gauss-Seidel迭代差分解图像') SOR迭代法function[u_3,m_3]=SOR_Solve(A,b,err,w,n)D=diag(A);D=diag(D);L=-tril(A,-1);R2=-triu(A,1);B1=(D-w*L)\((1-w)*D+w*R2);g=w*(D-w*L)\b;m_3=0;u_3_0=zeros(size(b,1),1);%初始迭代值flag=1;while flagu_3=B1*u_3_0+g;if norm(u_3-u_3_0,inf)<errflag=0;endu_3_0=u_3;m_3=m_3+1;%迭代次数值enduu=zeros(n+1);for mm=1:n-1for nn=1:n-1uu(nn+1,mm+1)=u_3((mm-1)*(n-1)+nn,1);endend%SOR迭代差分解图像x=[0:1/n:1];y=[0:1/n:1];figure(4)mesh(x,y,uu); title('SOR迭代差分解图像') PR方法clc;N=16;pi=3.1416;tao_opt=0.5*(1/sin(pi/N))%L1矩阵的定义L_1=zeros((N-1)*(N-1));for i=1:(N-1)*(N-1)-1L_1(i,i)=2;if mod(i,N-1)~=0L_1(i,i+1)=-1;L_1(i+1,i)=-1;endendL_1((N-1)^2,(N-1)^2)=2;%L2矩阵的定义L_2=zeros((N-1)*(N-1));for i=1:(N-1)*(N-1)-1L_2(i,i)=2;endL_2((N-1)^2,(N-1)^2)=2;for i=N:(N-1)^2L_2(i,i-(N-1))=-1;L_2(i-(N-1),i)=-1;end%右端项的定义b=zeros((N-1)^2,1);for j=1:N-1for i=1:N-1b((j-1)*(N-1)+i)=-(2*pi^2*exp(pi*(i/N+j/N))*sin(pi*(i/N+j/N)))/(N^2);endendI=eye((N-1)*(N-1));%PR格式的建立u=zeros((N-1)^2,1);for i=1:(N-1)^2u(i)=1;endu1=inv(I+tao_opt*L_1)*((I-tao_opt*L_2)*u+tao_opt*b);u2=inv(I+tao_opt*L_2)*((I-tao_opt*L_1)*u1+tao_opt*b);count=1;while norm(u2-u,inf)>=0.01count=count+1;u=u2;u1=inv(I+tao_opt*L_1)*((I-tao_opt*L_2)*u+tao_opt*b);u2=inv(I+tao_opt*L_2)*((I-tao_opt*L_1)*u1+tao_opt*b);endu2z=zeros((N-1)^2,1);for j=1:N-1for i=1:N-1z((j-1)*(N-1)+i)=exp(pi*(i/N+j/N))*(sin(pi*i/N)*sin(pi*j/N));endendzuu=zeros(N+1);for mm=1:N-1for nn=1:N-1uu(nn+1,mm+1)=u2((mm-1)*(N-1)+nn,1);endend%画图x=[0:1/N:1];y=[0:1/N:1];figure(1)mesh(x,y,uu);title('单参数PR迭代差分解图像（c1=c2=0.5*(1/sin(pi/16))') %画图数值解共轭梯度法function[u_5,m_5]=gongetidufa(A,b,n,err)m_5=0;u_5_0=zeros(size(b,1),1);p0=b-A*u_5_0;r0=p0;while flaga0=((r0)'*r0)/((p0)'*(A*p0));u_5=u_5_0+a0*p0;if norm(u_5-u_5_0,inf)<errflag=0;endu_5_0=u_5;r1=r0-a0*A*p0;b0=((r1)'*r1)/((r0)'*r0);p0=r1+b0*p0;r0=r1;m_5=m_5+1;enduu=zeros(n+1);for mm=1:n-1for nn=1:n-1uu(nn+1,mm+1)=u_5((mm-1)*(n-1)+nn,1);endend%共轭梯度法迭代差分解图像x=[0:1/n:1];y=[0:1/n:1];figure(5)mesh(x,y,uu);title('共轭梯度法迭代差分解图像')预处理共轭梯度法function[u_6,m_6]=yuchuligongetidufa(A,b,n,err) m_6=0;u_6_0=zeros(size(b,1),1);% D=diag(A);% D=diag(D);% L=tril(A,-1);% L=L+D\D;% B=L*D*(L)';B=0.7*eye((n-1)^2)p0=b-A*u_6_0;r0=p0;s0=B\r0;flag=1;while flaga0=((r0)'*s0)/((p0)'*(A*p0));u_6=u_6_0+a0*p0;if norm(u_6-u_6_0,inf)<errflag=0;endu_6_0=u_6;r1=r0-a0*A*p0;s1=B\r1;b0=((r1)'*s1)/((r0)'*s0);p0=s1+b0*p0;r0=r1;s0=s1;m_6=m_6+1;enduu=zeros(n+1);for mm=1:n-1for nn=1:n-1uu(nn+1,mm+1)=u_6((mm-1)*(n-1)+nn,1);endend%预处理共轭梯度法迭代差分解图像x=[0:1/n:1];y=[0:1/n:1];figure(6)mesh(x,y,uu);title('预处理共轭梯度法迭代差分解图像')五．实验数据分析1. 实验数据（部分）u精确u_1 u_2 u_3 u_4 u_5 u_6值0.05637 0.04884 0.04891 0.03396 0.04840 0.04901 0.049050.13455 0.11953 0.11967 0.08309 0.11869 0.11991 0.119910.23772 0.21531 0.21552 0.14965 0.21408 0.21585 0.215860.36820 0.33863 0.33891 0.23534 0.33708 0.33930 0.339350.52689 0.49057 0.49091 0.34090 0.48875 0.49139 0.491430.71247 0.67004 0.67043 0.46557 0.66802 0.67095 0.670940.92044 0.87287 0.87330 0.60645 0.87074 0.87383 0.873811.14208 1.09068 1.09112 0.75772 1.08851 1.09169 1.091631.36315 1.30963 1.31008 0.90978 1.30751 1.31062 1.310561.56264 1.50911 1.50954 1.04829 1.50711 1.51002 1.510001.71145 1.66037 1.66076 1.15331 1.65858 1.66118 1.661211.77124 1.72534 1.72568 1.19840 1.72384 1.72602 1.726071.69357 1.65570 1.65597 1.14999 1.65453 1.65624 1.656291.41964 1.39248 1.39267 0.96714 1.39166 1.39286 1.39289 0.88074 0.86649 0.86659 0.60181 0.86610 0.86670 0.86671 0.13455 0.11953 0.11967 0.08309 0.11869 0.11991 0.11991 0.32120 0.29120 0.29149 0.20241 0.28962 0.29194 0.29193 0.56747 0.52270 0.52313 0.36326 0.52036 0.52377 0.523780.87895 0.81982 0.82038 0.56969 0.81686 0.82120 0.821201.25777 1.18506 1.18574 0.82341 1.18158 1.18669 1.186681.70076 1.61573 1.61651 1.12257 1.61188 1.61752 1.617492.19723 2.10177 2.10261 1.46014 2.09769 2.10368 2.103612.72631 2.62302 2.62390 1.82216 2.61887 2.62499 2.624893.25404 3.14632 3.14721 2.18557 3.14227 3.14825 3.148173.73024 3.62234 3.62319 2.51612 3.61854 3.62413 3.624114.08547 3.98234 3.98313 2.76608 3.97894 3.98391 3.98397 4.22818 4.13538 4.13606 2.87229 4.13251 4.13672 4.13678 4.04279 3.96609 3.96663 2.75463 3.96386 3.96714 3.96720 3.38888 3.33378 3.33416 2.31541 3.33224 3.33453 3.33454 2.10245 2.07352 2.07372 1.44009 2.07277 2.07391 2.07391 0.23772 0.21531 0.21552 0.14965 0.21408 0.21585 0.215860.56747 0.52270 0.52313 0.36326 0.52036 0.52377 0.523781.00258 0.93567 0.93631 0.65019 0.93223 0.93721 0.937211.55288 1.46439 1.46523 1.01750 1.46004 1.46636 1.466392.22215 2.11316 2.11418 1.46817 2.10805 2.11550 2.115503.00480 2.87712 2.87828 1.99880 2.87146 2.87971 2.879693.88193 3.73830 3.73956 2.59692 3.73231 3.74106 3.741004.81668 4.66094 4.66226 3.23770 4.65485 4.66379 4.663685.74904 5.58628 5.58760 3.88030 5.58033 5.58908 5.589006.59037 6.42696 6.42823 4.46408 6.42138 6.42956 6.429547.21796 7.06143 7.06260 4.90462 7.05644 7.06374 7.063777.47011 7.32892 7.32993 5.09027 7.32471 7.33087 7.330937.14257 7.02561 7.02642 4.87949 7.02235 7.02718 7.027205.98727 5.90309 5.90365 4.09978 5.90082 5.90418 5.90417 2.Matlab作图H1 = 1/64时：精确解：Jacobi迭代解：Guass_Seidel迭代解：SOR迭代解：PR交替隐式差分迭代解：共轭梯度法解：预共轭梯度法解：3.各种迭代法性能分析h=1/64从上表可以看出，PR方法需要迭代次数最多，说明PR方法的收敛速度比较慢，预处理共轭梯度法需要迭代的次数最少，收敛较快；SOR方法虽然迭代的次数不多，但花费的时间却最长，说明该方法的效率最低。