实验名称微分方程数值解剖析
- 格式:doc
- 大小:351.50 KB
- 文档页数:11
常微分方程数值解实验报告学院:数学与信息科学专业:信息与计算科学姓名:郑思义学号: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尝试不同的步长,观察相应的数值解的变化。
数学中的微分方程数值解研究数学中的微分方程是一种描述自然界中变化的关系的数学工具,广泛应用于物理、工程、生物等领域。
微分方程数值解研究旨在通过使用计算机数值方法,求得微分方程的近似解,解决实际问题。
本文将介绍微分方程数值解的基本原理和常见方法。
一、微分方程数值解的原理微分方程是描述变化过程的数学模型,包含导数或微分项。
求解微分方程的精确解并不易于实现,因此需要借助数值方法来近似求解。
基本思路是将微分方程转化为离散的差分方程,通过计算机进行迭代计算,得到数值解。
二、常见的微分方程数值解方法1. 欧拉法欧拉法是最简单的数值解法之一,将微分方程转化为差分方程。
通过将自变量上的区间进行离散化,得到相应的插值点,然后利用差分逼近导数,进而求解差分方程。
欧拉法的优点是简单易懂,但对于复杂的微分方程可能会导致较大误差。
2. 改进的欧拉法改进的欧拉法是对欧拉法的改进,通过使用自变量在区间上的中点进行近似,减小误差。
改进的欧拉法是一种较为稳定和精确的数值解法,但对于某些特殊的微分方程,仍可能存在一定误差。
3. 龙格-库塔法龙格-库塔法是一种常用的数值解法,通过利用几个不同阶数的差分逼近,得到更精确的解。
常用的龙格-库塔法有四阶和二阶方法,其中四阶方法较为常用。
龙格-库塔法相较于前两种方法,具有更高的精确性和稳定性。
4. 有限差分法有限差分法是一种将微分方程转化为差分方程的数值解法。
通过将自变量上的区间分成若干个小区间,利用差商逼近导数,将微分方程转化为差分方程。
有限差分法适用于各种类型的微分方程,但需要选择合适的差分格式和网格,以保证数值解的准确性。
5. 贝尔斯托法贝尔斯托法是一种针对刚性微分方程的数值解法。
刚性微分方程是指在某些自变量的取值区间上,解的变化速度差别很大的微分方程。
贝尔斯托法通过选择合适的步长,以及使用隐式差分格式,对刚性微分方程进行数值求解,能够得到较为准确的结果。
三、微分方程数值解研究的应用微分方程数值解研究在实际问题中具有广泛的应用。
微分方程数值解法实验报告班级:姓名:学号:日期:一、实验目的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}误差-步长曲线如下:实验结论:通过对比欧拉方法和改进的欧拉方法的数值解与解析解的误差,可以发现改进的欧拉方法具有更高的精度和收敛性。
微分方程数值解实验报告微分方程数值解实验报告一、引言微分方程是数学中一类重要的方程,广泛应用于各个科学领域。
在实际问题中,往往难以得到微分方程的解析解,因此需要借助数值方法来求解。
本实验旨在通过数值解法,探索微分方程的数值解及其应用。
二、数值解法介绍常用的微分方程数值解法有欧拉法、改进欧拉法、四阶龙格-库塔法等。
在本实验中,我们将采用改进欧拉法进行数值解的求取。
改进欧拉法是一种一阶的显式迭代法,其基本思想是将微分方程的导数用差商来近似表示,并通过迭代逼近真实解。
具体迭代公式如下:\[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}\)。
将这些数值解进行可视化展示,可以更直观地观察到解的变化趋势。
根据实验结果,我们可以发现随着迭代次数的增加,数值解逐渐逼近了真实解。
探索实验8 常微分方程初值问题数值解一、 实验目的了解常微分方程初值问题的数值解概念,掌握解常微分方程初值问题的Euler 方法及改进的Euler 方法和Runge-Kutta 方法解常微分方程初值问题的算法构造和计算。
能用程序实现解常微分方程初值问题的Euler 方法、改进的Euler 方法和经典的Runge-Kutta 方法,学习用计算机求常微分方程初值问题数值解的一些科学计算方法和简单的编程技术。
二、概念与结论1. 常微分方程初值问题:常微分方程特解问题称为初值问题,通常其形式为2.常微分方程初值问题数值解:常微分方程初值问题的解y(x)在[a,b]上的有限个值y(x k )的近似值y k 称为常微分方程初值问题数值解,其中 x k = x k-1 + h k-1 ,k=1,2,3,…,N 。
x k 称为节点, h k 称为步长。
通常,步长h 不变,取为等距步长 h k =h=(b-a )/N ,N 为等分区间[a,b]分割数。
3.常微分方程初值问题数值解法:求常微分方程初值问题数值解y k 的方法称为微分方程初值问题数值解法。
4.单步法:在计算y k+1之时只用到y k 的方法,其计算公式有:显式单步法计算公式: y k+1=y k +h Φ(x k ,y k ;h )隐式单步法计算公式: y k+1 = y k +h Φ(x k ,y k ,y k+1,h )式中函数Φ是连续函数,称为增量函数。
5.多步法:在计算y k+1之时不仅用到y k ,还要用y k-1,y k-2,…,一般m 步法要用到y k ,y k-1,y k-2,… y k-m+1, 多步法也有显式方法和隐式方法之分。
6.数值解法的局部截断误差:假设某常微分方程初值问题数值解法在x= x k 没有误差,即y(x k )= y k ,称T k+1 =y(x k+1)- y k+1为该数值方法的局部截断误差。
显式单步法有其局部截断误差为:T k+1 =y(x k+1)- y(x k )- h Φ(x k ,y(x k ),h )7.数值解法的阶: 如果某常微分方程初值问题数值解法的局部截断误差为:T k+1 =O (h p+1)则称该数值方法的阶为P,或称该数值方法为P 阶方法。
数值方法的阶越高,方法越好。
三、程序中Mathematica 语句解释:1.k++,++k,k--,--k⎩⎨⎧=<<=')1()(),(0y a y b x a y x f yk++ 表示赋值关系k = k+1 ,如: k=1;Table[++k,{5}]获得表{2,3,4,5,6}++k 表示先处理k的值,再做赋值 k=k+1, 如: k=1;Table[k++,{5}]获得表{1,2,3,4,5} k-- 表示赋值关系 k = k-1, 如: k=1;Table[k--,{5}]获得表{1, 0, -1, -2, -3}--k 表示先处理k的值,再做赋值 k=k-1,如:k=1;Table[--k,{4}]获得表{0,-1,-2,-3} 2. x+=k, x*=kx+=k 表示 x = x + kx*=k 表示x = x * k3.For[stat,test,incr,body]以stat为初值,重复计算incr和body直到test为False终止。
这里start为初始值,test为条件,incr为循环变量修正式,body为循环体,通常由incr项控制test的变化。
注意: 上述命令形式中的start可以是由复合表达式提供的多个初值,如果循环体生成Break[ ] 语句,则退出For循环; 如果循环体生成Continue[ ] 语句,则由incr的增量进入For语句的下一次循环。
四、方法与程序在实际问题中,常常会遇到一些常微分方程初值问题。
对这些问题如果只求助于高等数学中介绍的用求通解加确定边界条件的方法去求解往往是无能为力的。
因为一方面通解可能求不出来,另一方面所求的解可能是不能用初等函数表达的形式。
人们处理这类问题主要采用常微分方程初值问题数值解的方法来求解。
这类方法有很多,这里主要介绍解常微分方程初值问题的Euler方法、改进的Euler方法和Runge-Kutta方法的构造过程和算法程序。
1.Euler方法Euler方法是最简单的求微分方程数值解的方法。
这个方法由于精度不高,实用中较少使用。
人们常用Euler方法来说明求微分方程数值解涉及到的一些问题。
1.1Euler方法的构造过程:Euler方法涉及的计算公式有很多方法,这里介绍在求微分方程数值解用的较多的Taylor展开法。
设f(x,y)充分光滑, x n+1=x n+h,将y(x n+1)在x n点作Taylor展开,得:y (x n+1)=y(x n)+hy'(x n)+(h2/2!)y"(x n)+O(h3)取其关于h 的线性部分,有:y (x n+1)≈y(x n)+hy'(x n)注意到y'(x n)= f(x n,y(x n)),用y k代替y(x k)并将“≈”换为等号“=”,得到Euler公式:y n+1= y n+h f(x n,y n)于是,由初始条件y0=y(a),借助Euler公式就可以依次计算出微分方程初值问题的数值解:y1 ,y2,…,y k称由Euler公式求数值解的方法为Euler方法。
显然Euler方法是单步显式方法。
因为对Euler公式有局部截断误差为T n+1= O(h2)因此Euler方法是一阶方法。
1.2 Euler方法算法:1.输入函数f(x,y)、初值y0、自变量区间端点a,b步长h2. 计算节点数n 和节点x k3. 用Euler 公式y n+1= y n +h f(x n ,y n ) 求数值解1.3 Euler 方法程序:Clear[x,y,f]f[x_,y_]= Input["函数f(x,y)="]y0=Input["初值y 0 ="]a=Input["左端点a="]b=Input["右端点b="]h=Input["步长h="]n=(b-a)/h;For[i=0,i<n,i++,xk=a+i*h;y1=y0+h*f[xk,y0];Print["y(",xk+h//N,")=",y1//N];y0=y1]说明:本程序用Euler 公式求常微方程初值问题数值解。
程序执行后,按要求通过键盘依次输入输入函数f(x,y)、初值y 0、自变量区间端点a,b 步长h 后,计算机则给出常微方程初值问题数值解。
程序中变量说明:f[x,y]: 存放函数f(x,y)y0: 存放初值y 0及数值解a :存放自变量区间左端点b: 存放自变量区间右端点n: 存放节点个数h: 存放节点步长xk :存放节点xiy1: 存放数值解注:1)语句Print["y(",xk+h//N,")=",y1//N]是将数值解用6位有效数字显示出来,如果要显示n 位数有效数字,可以将语句改为Print["y(",xk+h//N,")=",N[y1,n]]2)Mathematica 中有求微分方程初值问题数值解的命令,形式为:NDSolve[ {y’[x]==f[x,y[x]],y[a]==y0}, y, {x, a, b}]由N DSolve 命令得到的解是以{{未知函数名->InterpolatingFunction[range, <>]}}的形式给出的,其中的InterpolatingFunction[range, <>]是所求的插值函数表示的数值解, range 就是所求数值解的自变量范围。
1.4 例题与实验例1. 用Euler 方法求初值问题的数值解。
取步长h=0.1,并在一个坐标系中画出数值解与准确解的图形。
⎪⎩⎪⎨⎧=<<-='1)0(102y x y x y y解:执行Euler 方法程序后,在输入的窗口中按提示分别输入x-2x/y 、1、0、1、0.1,每次输入后用鼠标点击窗口的“OK ”按扭,计算机在屏幕上给出如下数值解结果:y(0.1)=1.1y(0.2)=1.19182y(0.3)=1.27744y(0.4)=1.35821y(0.5)=1.43513y(0.6)=1.50897y(0.7)=1.58034y(0.8)=1.64978y(0.9)=1.71778y(1.)=1.78477本题的准确解为y(x)= (1 + 2 x)1/2,在一个坐标系中画出数值解与准确解的图形如下:图中点图是数值解,曲线为准确解。
从图中可以知道数值解与准确解比较接近,但还是有误差的。
2. 改进的Euler 方法改进的Euler 方法比Euler 方法精度高,其中把微分方程初值问题数值隐式解法计算公式变为预测、校正公式的做法很有代表性。
2.1改进的Euler 方法构造过程:将微分方程y '=f(x,y) 在[x k ,x k+1]积分,得对右端的定积分用数值积分的梯形公式,有用y k 代替y(x k )并将“≈”换为等号“=”可得求微分方程初值问题的数值解的梯形公式:⎰+=-+11k k x x k k dxx y x f x y x y ))(,()()()))(,())(,((2))(,()()(1111++++≈=-⎰+k k k k x x k k x y x f x y x f h dx x y x f x y x y k k这个公式是单步隐式公式。
可以证明它是二阶方法,比Euler 方法阶高。
不过,在给定初始条件y 0=y(a)后,要求出数值解,每一次计算y k 的值都要进行非线性方程求根的迭代解法来完成,因此计算量很大。
为减少计算量,通常采用迭代一两次就转入下一步计算的方法,特别,如果先用Euler 公式计算一次以对要计算的下一步值解进行预测,然后再用梯形公式对其进行校正的方法得到下一步的值,它的计算格式为:这个计算格式称为改进的Euler 公式,而称由如上计算格式求数值解的方法为改进的Euler 方法。
2.2 改进的Euler 方法算法:1.输入函数f(x,y)、初值y 0、自变量区间端点a,b 步长h2.计算节点数n 和节点x k3.用改进的Euler 公式求数值解2.3 改进的Euler 方法程序:Clear[x,y,f]f[x_,y_]= Input["函数f(x,y)="]y0=Input["初值y 0 ="]a=Input["左端点a="]b=Input["右端点b="]h=Input["步长h="]n=(b-a)/h;For[i=0,i<n,i++,xk=a+i*h;y1=y0+h*f[xk,y0];y1=y0+h/2*(f[xk,y0]+f[xk+h,y1]);Print["y(",xk+h//N,")=",y1//N];y0=y1]注:语句h=Input["步长h="]的步长输入应该用带小数点的数输入,这样可以加快计算速度,否则,由于Mathematica 软件本身的原因,可能会出现计算时间过长等计算不出结果的情况。