常微分方程的求解与定性分析实验报告

常微分方程数值解实验报告学院：数学与信息科学专业：信息与计算科学姓名：郑思义学号：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进行计算，并与真解作比较。

常微分方程的数值解法专业班级：信息软件 姓名：吴中原 学号：120108010002 一、实验目的1、熟悉各种初值问题的算法，编出算法程序；2、明确各种算法的精度与所选步长有密切关系；通过计算更加了解各种 算法的优越性。

二、实验题目1、根据初值问题数值算法，分别选择二个初值问题编程计算；2、试分别取不同步长，考察某节点j x处数值解的误差变化情况； 3、试用不同算法求解某初值问题，结果有何异常； 4、分析各个算法的优缺点。

三、实验原理与理论基础（一） 欧拉法算法设计对常微分方程初始问题(6-1)(6-2)用数值方法求解时，我们总是认为(6-1)、(6-2)的解存在且唯一。

欧拉法是解初值问题的最简单的数值方法。

从(6-2)式由于y (x 0) = y 0已给定，因而可以算出),()('000y x f x y =。

设x 1 = h 充分小，则近似地有：),()(')()(00001y x f x y hx y x y =≈-(6-3)记 ,n ,,i x y y i i 10 )(== 从而我们可以取),(0001y x hf y y ==作为)(1x y 的近似值。

利用1y 及f (x 1, y 1)又可以算出)(2x y 的近似值：),(1112y x hf y y +=一般地，在任意点()h n x n 11+=+处)(x y 的近似值由下式给出),(1n n n n y x hf y y +=+(6-4)这就是欧拉法的计算公式，h 称为步长。

⎪⎩⎪⎨⎧==)( ),(d d 00y x y y x f x y（二）四阶龙格-库塔法算法设计：欧拉公式可以改写为：()111,i i i i y y k k hf x y +=+⎧⎪⎨=⎪⎩，它每一步计算(),f x y 的值一次，截断误差为()2o h 。

改进的欧拉公式可以改写为：()()()11212112,,i i i i i i y y k k k hf x y k hf x h y k +⎧=++⎪⎪=⎨⎪=++⎪⎩，它每一步要计算(),f x y 的值两次，截断误差为()3o h 。

实验报告实验项目名称常微分方程的数值解法实验室数学实验室所属课程名称微分方程数值解实验类型上机实验实验日期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 法的精度更高。

常微分方程实验报告一、实验目的常微分方程是数学分析和实际应用中非常重要的一部分，本次实验的主要目的是通过实际操作和计算，深入理解常微分方程的概念、性质和求解方法，并能够将其应用到实际问题中，提高我们解决数学问题和实际应用问题的能力。

二、实验原理常微分方程是指含有一个自变量和一个未知函数及其导数的等式。

求解常微分方程的方法有很多，常见的有变量分离法、一阶线性方程的求解方法（如常数变易法）、恰当方程的求解方法（通过积分因子）等。

对于一阶常微分方程，形如＼（y' ＋ p(x)y ＝ q(x)＼）的方程，可以使用积分因子＼（e^｛＼int p(x)dx}＼）来求解。

对于可分离变量的方程，形如＼（g(y)dy ＝ f(x)dx\），可以通过分别积分求解。

三、实验内容（一）一阶常微分方程的求解1、求解方程＼（y' ＋ 2xy ＝ 2x\）首先，计算积分因子＼（e^｛＼int 2xdx} ＝ e^｛x^2}＼），然后将方程两边乘以积分因子得到：＼（（ye^｛x^2}）＇＝ 2xe^｛x^2}＼）两边积分可得＼（ye^｛x^2} ＝ e^｛x^2} ＋ C\），解得＼（y ＝1 ＋ Ce^｛x^2}＼）2、求解方程＼（xy' y ＝ x^2\）将方程化为＼（y' ＼frac{y}｛x} ＝ x\），这里＼（p(x) ＝＼frac{1}｛x}＼），积分因子为＼（e^｛＼int ＼frac{1}｛x}dx} ＝＼frac{1}｛x}＼）。

方程两边乘以积分因子得到＼（（＼frac{y}｛x}）＇＝ 1\），积分可得＼（＼frac{y}｛x} ＝ x ＋ C\），即＼（y ＝ x^2 ＋ Cx\）（二）二阶常微分方程的求解1、求解方程＼（y'＇ 2y' ＋ y ＝ 0\）特征方程为＼（r^2 2r ＋ 1 ＝ 0\），解得＼（r ＝ 1\）（二重根），所以通解为＼（y ＝（C_1 ＋ C_2x)e^x\）2、求解方程＼（y'＇＋ 4y ＝ 0\）特征方程为＼（r^2 ＋ 4 ＝ 0\），解得＼（r ＝＼pm 2i\），所以通解为＼（y ＝ C_1\cos(2x) ＋ C_2\sin(2x)＼）（三）应用常微分方程解决实际问题1、考虑一个物体在受到与速度成正比的阻力作用下的运动，其运动方程为＼（m\frac{dv}｛dt} ＝ kv\）（其中＼（m\）为物体质量，＼（k\）为阻力系数），求解速度＼（v\）随时间＼（t\）的变化。

常微分方程数值解实验报告实验报告：常微分方程数值解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尝试不同的步长，观察相应的数值解的变化。

常微分方程实验报告《常微分方程》综合性实验实验报告实验班级05应数（3）学生姓名江晓荣学生学号200530770314指导老师方平华南农业大学理学院应用数学系实验微分方程在数学建模中的应用及数值解的求法一、实验目的1．了解常微分方程的基本概念。

2．常微分方程的解了解析解和数值解。

3．学习、掌握MA TLAB 软件有关求解常微分方程的解析解和数值解的有关命令。

4. 掌握微分方程在数学建模中的应用。

二、实验内容1．用MA TLAB 函数dsolve 符号求解常微分方程的通解和特解。

2．用MA TLAB 软件数值求解常微分方程。

三、实验准备1．用MA TLAB 求常微分方程的解析解的命令用MA TLAB 函数dsolve 求常微分方程()(,,,,,,)0n F x y y y y y ''''''= （7.1）的通解的主要调用格式如下：S=dsolve('eqn', 'var')其中输入的量eqn 是改用符号方程表示的常微分方程(,,,2,)0F x y Dy D y Dny = ，导数用D 表示，2阶导数用D2表示，以此类推。

var 表示自变量，默认的自变量为t 。

输出量S 是常微分方程的解析通解。

如果给定常微分方程（7.1）的初始条件()00010(),(),,()n n y x a y x a y x a '=== ，则求方程（7.1）的特解的主要调用格式如下：S=dsolve('eqn', 'condition1 ',…'conditonn ','var')其中输入量eqn ，var 的含义如上，condition1,…conditonn 是初始条件。

输出量S 是常微分方程的特解。

2．常微分方程的数值解法除常系数线性微分方程可用特征根法求解、少数特殊方程可用初等积分法求解外，大部分微分方程无解析解，应用中主要依靠数值解法。

1. 分别用Euler 法和ode45解下列常微分方程并与解析解比较： (1) ,(0)1,13y x y y x '=+=<<Euler 法：function [t,y]=euler(Fun,tspan,y0,h) t=tspan(1):h:tspan(2); y(1)=y0;for i=1:length(t)-1y(i+1)=y(i)+h.*feval(Fun,t(i),y(i)); end t=t'; y=y';function f=Fun(x,y) % 常微分方程的右端函数 f=x+y;>> [x,y]=euler('Fun',[0,3],1,0.1)>> [x,y] ans =0 1.0000 0.1000 1.1000 0.2000 1.2200 0.3000 1.3620 0.4000 1.5282 0.5000 1.7210 0.6000 1.9431 0.7000 2.1974 0.8000 2.4872 0.9000 2.8159 1.0000 3.1875 1.1000 3.6062 1.2000 4.0769 1.3000 4.6045 1.4000 5.1950 1.5000 5.8545 1.6000 6.5899 1.7000 7.4089 1.8000 8.3198 1.9000 9.3318 2.0000 10.4550 2.1000 11.7005 2.2000 13.0805 2.3000 14.6086 2.4000 16.2995 2.5000 18.1694 2.6000 20.2364 2.7000 22.5200 2.8000 25.0420 2.9000 27.8262 3.0000 30.8988ode45：>> [x,y]=ode45('Fun',[0,3],1) ans =0 1.0000 0.0502 1.0528 0.1005 1.1109 0.1507 1.17460.2010 1.2442 0.2760 1.3596 0.3510 1.4899 0.4260 1.63610.5010 1.7996 0.5760 1.9817 0.6510 2.1838 0.7260 2.40740.8010 2.6544 0.8760 2.9264 0.9510 3.2254 1.0260 3.55351.1010 3.9131 1.1760 4.3065 1.2510 4.7364 1.3260 5.20561.4010 5.7172 1.4760 6.2744 1.5510 6.8810 1.6260 7.54061.7010 8.2574 1.7760 9.0359 1.8510 9.8808 1.9260 10.79742.0010 11.7912 2.0760 12.8683 2.1510 14.0351 2.2260 15.29862.3010 16.6664 2.3760 18.1466 2.4510 19.7478 2.5260 21.47962.6010 23.3522 2.6760 25.3764 2.7510 27.5641 2.8260 29.92812.9010 32.4820 2.9257 33.3694 2.9505 34.2796 2.9752 35.21343.0000 36.1711解析解：>> y=dsolve('Dy=x+y','y(0)=1','x') y =2*exp(x) - x - 1(2) 20.01()2sin(),(0)0,(0)1,05y y y t y y t ''''-+===<< Euler 法：function f=Fun(t,y)% 常微分方程的右端函数f=[y(2);0.01*y(2)^2-2*y(1)+sin(t)];>> [t,y]=euler('Fun',[0,5],[0,1],0.2) ode45：>> [t,y]=ode45('Fun',[0,5],[0,1])t =0 0.0001 0.0001 0.0002 0.0002 0.0005 0.0007 0.0010 0.0012 0.00250.0037 0.0050 0.0062 0.0125 0.0188 0.0251 0.0313 0.0627 0.0941 0.12550.1569 0.2819 0.4069 0.5319 0.6569 0.7819 0.9069 1.0319 1.1569 1.28191.4069 1.5319 1.6569 1.7819 1.90692.0319 2.1569 2.2819 2.4069 2.53192.6569 2.7819 2.90693.0319 3.1569 3.2819 3.4069 3.5319 3.6569 3.78193.90694.0319 4.1569 4.2819 4.4069 4.5319 4.6569 4.7427 4.8285 4.91425.0000y =0 1.0000 0.0001 1.0000 0.0001 1.0000 0.0002 1.0000 0.0002 1.00000.0005 1.0000 0.0007 1.0000 0.0010 1.0000 0.0012 1.0000 0.0025 1.00000.0037 1.0000 0.0050 1.0000 0.0062 1.0000 0.0125 1.0000 0.0188 1.00000.0251 0.9999 0.0313 0.9998 0.0627 0.9987 0.0941 0.9965 0.1253 0.99340.1564 0.9893 0.2786 0.9632 0.3966 0.9220 0.5085 0.8662 0.6126 0.79670.7072 0.7146 0.7908 0.6210 0.8620 0.5176 0.9198 0.4058 0.9632 0.28760.9915 0.1647 1.0043 0.0392 1.0013 -0.0869 0.9826 -0.2117 0.9485 -0.33310.8996 -0.4490 0.8365 -0.5578 0.7605 -0.6577 0.6725 -0.7471 0.5742 -0.8246 0.4669 -0.8889 0.3525 -0.9393 0.2327 -0.9748 0.1095 -0.9950 -0.0154 -0.9996 -0.1398 -0.9887 -0.2619 -0.9624 -0.3798 -0.9212 -0.4916 -0.8657 -0.5957 -0.7970 -0.6904 -0.7161 -0.7742 -0.6242 -0.8460 -0.5228 -0.9046 -0.4134 -0.9491 -0.2978 -0.9789 -0.1777 -0.9934 -0.0549 -0.9945 0.0300 -0.9883 0.1146 -0.9748 0.1985 -0.9543 0.28092. 求一通过原点的曲线，它在(,)x y 处的切线斜率等于22,0 1.57.x y x +<<若x 上限增为1.58，1.60会发生什么？function f=Fun(x,y) % 常微分方程的右端函数 f=2*x+y.^2;>> [x,y]=ode45('Fun',[0,1.57],0) x =0 0.0393 0.0785 0.1178 0.1570 0.1963 0.2355 0.2748 0.3140 0.3533 0.3925 0.4318 0.4710 0.5103 0.5495 0.5888 0.6280 0.6673 0.7065 0.7458 0.7850 0.8243 0.8635 0.9028 0.9420 0.9813 1.0205 1.0598 1.0990 1.1383 1.1775 1.2168 1.2560 1.2953 1.3345 1.3738 1.4130 1.4248 1.4367 1.4485 1.4604 1.4722 1.4840 1.4959 1.5077 1.5140 1.5203 1.5265 1.5328 1.5376 1.5424 1.5472 1.5519 1.5543 1.5567 1.5591 1.5614 1.5631 1.5647 1.5664 1.5681 1.5685 1.5690 1.5695 1.5700 y =0 0.0015 0.0062 0.0139 0.0247 0.0386 0.0556 0.0758 0.09920.1259 0.1559 0.1895 0.2266 0.2675 0.3124 0.3615 0.4152 0.4738 0.5378 0.6076 0.6841 0.7679 0.8601 0.9620 1.0751 1.2014 1.3434 1.5045 1.6892 1.9037 2.1557 2.4577 2.8282 3.3003 3.9056 4.7317 5.9549 6.4431 7.0116 7.6832 8.4902 9.4821 10.7170 12.3090 14.4551 15.9220 17.7080 19.9390 22.8164 25.6450 29.2282 33.9673 40.5910 44.9434 50.3088 57.1229 66.1087 74.3108 84.7123 98.4901 117.7875 124.9206 132.9699 142.1268 152.6415若x 上限增为1.58，1.60,则超出运算的范围，发生溢出。

常微分方程的定性分析常微分方程是研究自变量只涉及一个变量的微分方程，在科学和工程中具有广泛的应用。

定性分析是常微分方程中重要的一部分，它是指通过分析方程的性质和图像来揭示方程的解的行为。

在本文中，我们将讨论常微分方程的定性分析的基本方法和技巧。

一、平衡点和稳定性分析在进行定性分析之前，首先需要确定方程的平衡点。

平衡点是指微分方程中导数为零的点，即解保持恒定的点。

通过求解方程等于零的情况，我们可以找到方程的平衡点。

确定平衡点后，我们需要分析平衡点的稳定性。

稳定性是指当初始条件接近平衡点时，解是否会趋向于平衡点。

通过线性化的方法可以分析平衡点的稳定性，即在平衡点附近做泰勒展开，然后分析展开式的特征根。

二、相图和相轨线相图是用来描述微分方程解的整体行为的图形表示。

在相图中，自变量通常表示时间，因变量表示微分方程的解。

通过绘制相图，我们可以看到解的轨迹和相位变化。

相轨线是相图中的曲线，表示微分方程解在相空间中的轨迹。

通过绘制相轨线，我们可以直观地了解方程的解的行为。

相轨线可以通过数值方法或者解析方法进行求解。

三、参数分析和稳定性改变在定性分析中，我们可以通过改变微分方程中的参数来观察解的行为的变化。

通过参数的分析，我们可以看到解在不同参数取值下的定性变化。

特别是可以通过稳定性分析，观察参数的改变对平衡点的稳定性有何影响。

四、存在性和唯一性在进行定性分析之前，我们需要先讨论微分方程解的存在性和唯一性问题。

存在性指的是在给定的初始条件下是否存在解。

唯一性指的是解是否是唯一的。

通过利用积分器的理论可以证明微分方程解的存在性和唯一性。

五、应用实例下面通过几个实例来说明常微分方程定性分析的具体应用。

例1：考虑简谐振动方程m*x''+c*x'+k*x=0。

分析方程的解的稳定性和相轨线。

解：首先确定平衡点。

当加速度为零时，m*x''+c*x'+k*x=0，可得平衡点为x=0。

学生实验报告（4）一、实验综述. 归纳和学习求解常微分方程(组)的基本原理和方法；2. 掌握解析、数值解法，并学会用图形观察解的形态和进行解的定性分析；3. 熟悉MATLAB软件关于微分方程求解的各种命令；4. 通过范例学习建立微分方程方面的数学模型以及求解全过程；通过该实验的学习，使学生掌握微分方程(组)求解方法（解析法、欧拉法、梯度法、改进欧拉法等），对常微分方程的数值解法有一个初步了解，同时学会使用MATLAB软件求解微分方程的基本命令，学会建立微分方程方面的数学模型。

这对于学生深入理解微分、积分的数学概念，掌握数学的分析思维方法，熟悉处理大量的工程计算问题的方法是十分必要的。

二、实验过程（实验步骤、记录、数据、分析）1．开启MATLAB软件平台，开启MATLAB编辑窗口；2．根据问题，建立的线性规划模型，并编写求解规划模型的M文件；3．保存文件并运行；4．观察运行结果(数值或图形)，并不断地改变参数设置观察运行结果；5．根据观察到的结果和体会，写出实验报告。

三、实验要求与任务根据实验内容和步骤，完成以下实验，要求写出实验报告（实验目的→问题→数学模型→算法与编程→计算结果→分析、检验和结论）1．求微分方程的解析解，并画出它们的图形。

y '= y + 2 x, y (0) = 1, 0< x <1；程序如下：由y=dsolve('Dy=y+2*x','y(0)=1','x')得出解析解y =-2*x-2+3*exp(x) 建立函数m 文件：function y=myfun4(x)y=-2*x-2+3*exp(x)画图函数为fplot('myfun4',[0,1])图形如下：2．求微分方程⎪⎩⎪⎨⎧====-+]100[0)0(;0)0(01.03t uu u u u 的数值解，要求编写求解程序。

首先建立函数M 文件：function dy=myfun5(t,y)dy=zeros(2,1);dy(1)=y(2);dy(2)=0.1*(y(1).^3)-y(1);输入命令：[T,Y]=ode15s('myfun5',[0,10],[0,1]);plot(T,Y(:,1),'-',T,Y(:,2),'*');图形如下：3．Rossler 微分方程组：⎪⎩⎪⎨⎧-+=+=--=)('''c x z b z ayx y z y x 当固定参数b =2，c =4时，试讨论随参数a 由小到大变化（如 a ∈(0，0.65))而方程解的变化情况，并且画出空间曲线图形，观察空间曲线是否形成混沌状？ 首先建立函数M 文件rossler.m ，在其中用x(1)表示x ，用x(2)表示y ，用x(3)表示z.function r=rossler(t,x)global a;global b;global c;r=[-x(2)-x(3);x(1)+a*x(2);b+x(3)*(x(1)-c)];主程序如下：global a;global b;global c;b=2;c=4;t0=[0,200];for a=0:0.02:0.65[t,x]=ode45('rossler',t0,[0,0,0]);asubplot(1,2,1);plot(t,x(:,1),'r',t,x(:,2),'g',t,x(:,3),'b');title('x(红色),y(绿色),z(篮色)随t变化情况');xlabel('t'); subplot(1,2,2);plot3(x(:,1),x(:,2),x(:,3))title('相图');xlabel('x');ylabel('y');zlabel('z');pauseend当a=0 时，(x,y,z)收敛于（0,0.5,0.5）当a=0.05 时，(x,y,z)仍然收敛，但收敛速度较小。

常微分方程课程实验报告实验名称 Matlab在常微分方程中的应用z =C2*exp(-2*t)+C3*exp(2*t)【3．1】作图(先做出x关于t，y关于t,z关于t的函数图像)：hold on;for c1=0:0.1:1%C1、C2和C3都大于0for c2=0:0.1:1for c3=0:0.1:1t=0:0.1:1;x=c1.*exp(2.*t)+c3.*exp(-t);y=c1.*exp(2.*t)+c2.*exp(-t)+exp(-2.*t)*c3;z=c2.*exp(2.*t)+exp(-2.*t)*c3;subplot(1,3,1),plot(t,x),legend('t~x') %x关于t的函数图像 subplot(1,3,2),plot(t,y),legend('t~y') % y关于t的函数图像 subplot(1,3,3),plot(t,z),legend('t~z') % z关于t的函数图像endendend【3.2】作图（再做出x、y、z的三维图像）：hold on;X=[];Y=[];Z=[];for c1=0:0.2:1%C1、C2和C3都大于0for c2=0:0.2:1for c3=0:0.2:1t=0:0.1:1;x=c1.*exp(2.*t)+c3.*exp(-t);y=c1.*exp(2.*t)+c2.*exp(-t)+exp(-2.*t)*c3;z=c2.*exp(2.*t)+exp(-2.*t)*c3;X=[X,x];Y=[Y,y];Z=[Z,z];endColumns 37 through 410.9753 1.0129 1.0521 1.0929 1.1353第4题：【1】编写函数M文件——cwf1.m(用于求近似解)function dy=cwf1(x,y)dy=(cos(x)-2*x*y)/(x^2-1)【2】编写脚本M文件——chang8.m[X,Y]=ode45('cwf1',[0,1],1)y1=Y’ %转置【3】则求出的近似解为（x取值为0~1）：>>y1=Columns 1 through 101.0000 0.9756 0.9524 0.9303 0.9093 0.8892 0.8701 0.8520 0.8347 0.8183Columns 11 through 200.8028 0.7880 0.7742 0.7611 0.7488 0.7374 0.7269 0.7172 0.7085 0.7008Columns 21 through 300.6941 0.6886 0.6843 0.6815 0.6802 0.6809 0.6837 0.6890 0.6976 0.7102Columns 31 through 400.7277 0.7519 0.7851 0.8319 0.8964 0.9910 1.1404 NaN NaN NaNColumn 41-Inf【4】在command窗口运行以下语句，用于求精确解：>> y=dsolve('(x^2-1)*Dy+2*x*y-cos(x)=0','y(0)=1','x') %精确解。

微分方程实验报告实验名称：不同数值方法解常微分方程的数值精度实验目的：比较同一数值方法不同刨分数解常微分方程的数值精度,比较不同数值方法在相同刨分数下解常微分方程的数值精度,并将他们与理论数值精度相比；求出截断误差。

实验要求：1.编写程序求解微分方程000(,)()du f x u x x b dx u x u ⎧=⎪<≤⎨⎪=⎩其中f 是x 和u 的已知函数，0u 是给定初值。

2.至少用下列方法实验步骤：一. 初值问题：(),5.11,11,222≤<⎪⎩⎪⎨⎧=-='x u u xx u u 解析解：()3122334⎪⎪⎭⎫ ⎝⎛-=x x x u .二. 解：对该问题，10=u ，00=x ，mh x m =。

0125.0=h1. Euler 公式：1(,),0,1,m m m m u u hf x u m +=+=2. 梯形法公式：111[(,)(,)]/2,0,1,m m m m m m u u h f x u f x u m +++=++=3. 三阶Heun 公式：21131132(3)/4(,),(/3,/3),(2/3,2/3),0,1,m m m m m m m m u u hf k k k f x u k f x h u hk k f x h u hk m +=++==++=++=4. 四阶Runge-Kutta 公式：211234113243(22)/6(,),(/2,/2),(/2,/2),(,),0,1,m m m m m m m m m u u hf k k k k k f x u k f x h u hk k f x h u hk k f x h u hk m +=++++==++=++=++=5. 四阶Adams 外插公式：（显式）),(,4,3,24/)9375955/(3211m m m m m m m m m u x f f m f f f f h u u ==-+-+=---+6. 四阶Gear 公式：11231(4836163)/2512/25,3,4,m m m m m m u u u u u hf m +---+=-+-+=理论数值精度p=4；其中，m u 作为()m u x 的近似值（m=1,2,…）三．计算方法的数值精度 首先取步长h 进行计算，得到()u x 在b 点的近似值()h u b ，根据误差估计，得()()ph h e u b u b M h =-≈，其中p 是方法的阶。

常微分方程的数值解法实验报告实验报告：常微分方程的数值解法摘要：常微分方程（ODE）是描述动力学系统中物理量随时间变化的数学方程，广泛应用于自然科学和工程领域。

然而，对于一些复杂的非线性ODE，很难找到解析解。

因此，我们需要数值解法来求解这些方程。

本实验报告将介绍四种常见的常微分方程数值解法，分别是欧拉法、改进的欧拉法、四阶龙格-库塔法和自适应步长的龙格-库塔法，并通过数值实验比较它们的精度和效率。

1.引言在实际问题中，许多物理量的变化规律可以由常微分方程描述。

然而，对于复杂的非线性ODE，很难找到解析解。

因此，为了解决这类问题，我们需要借助数值方法来求解。

2.方法本实验采用四种常见的常微分方程数值解法：欧拉法、改进的欧拉法、四阶龙格-库塔法和自适应步长的龙格-库塔法。

（1）欧拉法是最简单的数值解法，通过将微分方程转化为差分方程，使用离散的步长来近似微分方程。

（2）改进的欧拉法在欧拉法的基础上进行了改进，使用预估-校正的方法来提高精度。

（3）四阶龙格-库塔法是一种经典的数值解法，通过利用不同步长处的斜率来近似微分方程，具有较高的精度。

（4）自适应步长的龙格-库塔法是在四阶龙格-库塔法的基础上改进而来的，根据步长的大小自适应地选择不同的步长，同时保证精度和效率。

3.实验设计为了比较这四种数值解法的精度和效率，我们设计了两个实验。

实验一是求解一阶常微分方程：dy/dx = -2x，初始条件y(0) = 1，解析解为y = 1 - x^2、实验二是求解二阶常微分方程：d^2y/dx^2 + y = 0，初始条件y(0) = 0，dy/dx(0) = 1，解析解为y = sin(x)。

4.结果与分析实验一中，比较四种数值解法在不同步长下的近似解和解析解，计算其误差。

实验结果表明，四阶龙格-库塔法和自适应步长的龙格-库塔法具有较高的精度，而欧拉法和改进的欧拉法的精度较低。

实验二中，我们比较四种数值解法在不同步长下的近似解和解析解，并计算其误差。

常微分方程实验报告(一)实验名称 姓名： 学 号 班级：一、 实验目的1. 学会用工具函数dsolve()求微分方程的解析解；2. 学会用ode45,ode23求微分方程的数值解。

二、实验原理1. 利用函数dsolve()计算常微分方程的符号解，格式为y=dsolve(‘方程1’,’方程2’,…,’初始条件1’,’初始条件2’,…,’自变量’) 自变量缺省值为t ，导数用D 表示，二阶导数用D2表示，以此类推，y 反回方程的解析解。

2. 利用ode45( ),ode23( )计算常微分方程的数值解，格式为[t,y]=ode45(‘odefun ’,[t0,tf],y0)采用变步长四节Runge-Kutta 法和五阶Runge-Kutta-Felhberg 法求数值解。

odefun:微分方程； t0：自变量的初值； tf ：自变量的终值； y0：初始向量值。

输出向量t 表示节点01(,,...,)n t t t 输出矩阵y 表示对应i t 的数值解。

三、 实验范例1. 求下列微分方程的解析解 (1) y ay b '=+y=dsolve(‘Dy=a*y+b ’,’x ’)(2) sin 2,(0)0,(0) 1.y y y x y y ''''++===s=dsolve(‘D2y+Dy+y=sin(2*x)’,’y(0)=0,Dy(0)=1’,’x’) (3) ,,(0)1,(0)1f f g g g f f g ''''=+=-==s=dsolve(‘Df=f+g,Dg=g-f ’,’Df(0)=1,Dg(0)=1’,’t ’)2. 求下列方程的数值解 (1) 1,(0) 1.y y x y '=-++=odefun=inline(‘-y+x+1’,’x ’,’y ’); [x,y]=ode45(odefun,[0,1],1)(2) 123213312()()()()0.51()()y y t y t y y t y t y y t y t '=⎧⎪'=-⎨⎪'=-⎩在[0,12]内的数值解,且满足初始条件123(0)0(0)1(0)1y y y =⎧⎪=⎨⎪=⎩。

常微分方程定性分析常微分方程（Ordinary Differential Equation, ODE）是数学领域的重要理论，用于描述变量及其导数之间的关系。

在实际应用中，常微分方程可以帮助我们理解自然现象、物理规律及经济现象等，因此对于常微分方程的定性分析具有重要意义。

一、常微分方程的基本概念和分类常微分方程是指未知函数的导数只涉及一个独立变量的方程。

一般形式为：dy/dx = f(x, y)其中，y是未知函数，x是自变量，f(x, y)是已知函数。

常微分方程可以分为一阶和高阶，线性和非线性等多种类型。

二、定性分析方法常微分方程的定性分析是指通过研究方程的性质和解的行为，确定方程解的大致形态和特征。

常用的定性分析方法有以下几种：1. 平衡点和稳定性分析平衡点是指满足dy/dx = 0的解点。

通过计算f(x, y)在平衡点处的导数，可以判断平衡点的稳定性。

若f'(x, y) > 0，则平衡点是不稳定的；若f'(x, y) < 0，则平衡点是稳定的。

2. 相图分析相图是指将导数关于自变量和因变量绘制成的图形。

根据相图的形态，可以初步判断方程解的行为。

常见的相图形态有稳定点、周期解、不稳定点等。

3. 变量分离法对于一些特殊形式的常微分方程，可以利用变量分离法进行求解。

变量分离法是指将方程中的自变量和因变量分开，再进行积分求解。

4. 拉普拉斯变换拉普拉斯变换是一种将常微分方程转化为代数方程的方法，通过对方程应用拉普拉斯变换，可以得到方程的解析解。

三、应用实例常微分方程定性分析在实际问题中具有广泛的应用。

以生态学模型为例，生态学家研究生物种群的增长与环境因素之间的关系时，常常采用常微分方程来描述种群数量的变化。

通过定性分析可以评估种群的稳定性，分析环境因素对种群数量的影响。

另外，常微分方程定性分析还可以应用于控制理论、金融学等领域。

在控制理论中，常微分方程用于描述动态系统的演化过程，通过定性分析可以预测系统的行为并设计相应的控制策略。