常微分方程上机实验报告
常微分方程 课程实习报告
福建农林大学计算机与信息学院（数学类课程）课程实习报告课程名称：常微分方程课程实习实习题目：常微分方程数值求解问题的实习姓名：上官火玉系：信息与计算科学专业：信息与计算科学年级：2008学号：081152048指导教师：陈永雪职称：讲师2010 年12 月 1 日福建农林大学计算机与信息学院数学类课程实习报告结果评定目录1. 实习的目的和任务 (1)2. 实习要求 (1)3. 实习地点 (1)4. 主要仪器设备 (1)5. 实习内容 (1)5.1 用不同格式对同一个初值问题的数值求解及其分析 (1)5.1.1求精确解 (1)5.1.2用欧拉法求解 (3)5.1.3用改进欧拉法求解 (5)5.1.4用4级4阶龙格—库塔法求解 (7)5.1.5 问题讨论与分析 (10)5.2 一个算法不同不长求解同一个初值问题及其分析 (12)6. 结束语 (13)参考文献 (13)常微分方程课程实习1.实习的目的和任务目的：通过课程实习能够应用MATLAB软来计算微分方程（组）的数值解；了解常微分方程数值解。
任务：通过具体的问题，利用MATLAB软件来计算问题的结果，分析问题的结论。
2.实习要求能够从案例的自然语言描述中，抽象出其中的数学模型；能够熟练应用所学的数值解计算方法；能够熟练使用MATLAB软件；对常微分方程数值解有所认识，包括对不同算法有所认识和对步长有所认识。
3.实习地点数学实验室4.主要仪器设备计算机、Microsoft Windows XPMatlab 6.55.实习内容5.1 用欧拉方法，改进欧拉方法，4阶龙格—库塔方法分别求下面微分方程的初值dy/dx=y+x+1 y(0)=1 x∈[0，2]5.1.1求精确解首先可以求得其精确解为：y=3*exp(x)-x-25.1.1 程序代码：>> x=0:0.1:2;>> y=3*exp(x)-x-2>> plot(x,y,'b*-');>> Data=[x',y']y =Columns 1 through 91.0000 1.2155 1.4642 1.74962.0755 2.44622.86643.3413 3.8766Columns 10 through 184.47885.1548 5.91256.76047.70798.76569.9451 11.2591 12.7218Columns 19 through 2114.3489 16.1577 18.1672Data =0 1.00000.1000 1.21550.2000 1.46420.3000 1.74960.4000 2.07550.5000 2.44620.6000 2.86640.7000 3.34130.8000 3.87660.9000 4.47881.0000 5.15481.1000 5.91251.2000 6.76041.3000 7.70791.4000 8.76561.5000 9.94511.6000 11.25911.7000 12.72181.8000 14.34891.9000 16.15772.0000 18.1672>>00.20.40.60.81 1.2 1.4 1.6 1.825.1.2 用欧拉法求解程序如下：建立函数文件cwfa1.mfunction [x,y]=cwfa1(fun,x_span,y0,h)x=x_span(1):h:x_span(2);y(1)=y0;for n=1:length(x)-1y(n+1)=y(n)+h*feval(fun,x(n),y(n)); endx=x';y=y';在MATLAB输入以下程序：>> clear all>> fun=inline(' y+x+1');>> [x,y]=cwfa1(fun,[0,2],1,0.1);>> [x,y]>> plot(x,y,'r*-')结果及其图象：ans =0 1.00000.1000 1.20000.2000 1.43000.3000 1.69300.4000 1.99230.5000 2.33150.6000 2.71470.7000 3.14620.8000 3.63080.9000 4.17381.0000 4.78121.1000 5.45941.2000 6.21531.3000 7.05681.4000 7.99251.5000 9.03171.6000 10.18491.7000 11.46341.8000 12.87981.9000 14.44772.0000 16.182500.20.40.60.81 1.2 1.4 1.6 1.825.1.3用改进欧拉法求解：程序如下：建立函数文件cwfa2.mfunction [x,y]=cwfa2(fun,x_span,y0,h)x=x_span(1):h:x_span(2);y(1)=y0;for n=1:length(x)-1k1=feval(fun,x(n),y(n));y(n+1)=y(n)+h*k1;k2=feval(fun,x(n+1),y(n+1));y(n+1)=y(n)+h*(k1+k2)/2;endx=x';y=y';在MATLAB输入以下程序：>> clear all>> fun=inline(' y+x+1');>> [x,y]=cwfa2(fun,[0,2],1,0.1); >> [x,y]>> plot(x,y,'r+-')结果及其图象：ans =0 1.00000.1000 1.21500.2000 1.46310.3000 1.74770.4000 2.07270.5000 2.44230.6000 2.86130.7000 3.33470.8000 3.86840.9000 4.46851.0000 5.14221.1000 5.89721.2000 6.74191.3000 7.68581.4000 8.73931.5000 9.91391.6000 11.22241.7000 12.67871.8000 14.29851.9000 16.09882.0000 18.0987>>00.20.40.60.81 1.2 1.4 1.6 1.825.1.4 用4阶龙格—库塔求解程序如下：建立函数文件cwfa3.mfunction [x,y]=cwfa3(fun,x_span,y0,h)x=x_span(1):h:x_span(2);y(1)=y0;for n=1:length(x)-1k1=feval(fun,x(n),y(n));k2=feval(fun,x(n)+h/2,y(n)+h/2*k1);k3=feval(fun,x(n)+h/2,y(n)+h/2*k2);k4=feval(fun,x(n+1),y(n)+h*k3);y(n+1)=y(n)+h*(k1+2*k2+2*k3+k4)/6;endx=x';y=y';在MATLAB输入以下程序：>> clear all;>> fun=inline(' y+x+1');>> [x,y]=cwfa3(fun,[0,2],1,0.1); >> [x,y]>> plot(x,y, 'b+-')结果及其图象：ans =0 1.00000.1000 1.21550.2000 1.46420.3000 1.74960.4000 2.07550.5000 2.44620.6000 2.86640.7000 3.34130.8000 3.87660.9000 4.47881.0000 5.15481.1000 5.91251.2000 6.76031.3000 7.70791.4000 8.76561.5000 9.94511.6000 11.25911.7000 12.72181.8000 14.34891.9000 16.15772.0000 18.1671>>00.20.40.60.81 1.2 1.4 1.6 1.825.1.5 问题讨论与分析由以上数值分析结果绘制表格：x=0:0.1:2;y1 = [1.0000 1.2155 1.4642 1.7496 2.0755 2.4462 2.8664 3.3413 3.8766 4.4788 5.1548 5.9125 6.7604 7.7079 8.7656 9.9451 11.2591 12.7218 14.3489 16.1577 18.1672];>> y1=[1.0000 1.2155 1.4642 1.7496 2.0755 2.4462 2.8664 3.3413 3.8766 4.4788 5.1548 5.9125 6.7604 7.7079 8.7656 9.9451 11.2591 12.7218 14.3489 16.1577 18.1672];>> y2=[1.0000 1.2000 1.4300 1.6930 1.9923 2.3315 2.7147 3.1462 3.6308 4.1738 4.7812 5.4594 6.2153 7.0568 7.9925 9.0317 10.1849 11.4634 12.8798 14.4477 16.1825];>> y3=[1.0000 1.2150 1.2192 1.4631 1.7477 2.0727 2.4423 2.8613 3.8684 4.4685 5.1422 5.8972 6.7419 7.6858 8.7393 9.9139 11.2224 12.6787 14.2985 16.0988 18.0987];>> y4=[1.0000 1.2155 1.4642 1.7496 2.0755 2.4462 2.8664 3.3413 3.8766 4.4788 5.1548 5.9125 6.7603 7.7079 8.7656 9.9451 11.2591 12.7218 14.3489 16.1577 18.1671 ];>> plot(x,y1,'r+-')>> hold on,plot(x,y2,'b-')>> plot(x,y1,'r+-')>> hold on,plot(x,y3,'b-')>> plot(x,y1,'r+-')>> hold on,plot(x,y4,'b-')00.20.40.60.81 1.2 1.4 1.6 1.8200.20.40.60.81 1.2 1.4 1.6 1.8200.20.40.60.81 1.2 1.4 1.6 1.82由上表和图可以看出欧拉法误差最大，而改进欧拉和龙格—库塔方法误差相对较小，并且龙格—库塔方法误差最小且大部分值都跟精确值相同。
实验报告——常微分方程的数值解法
实验报告实验项目名称常微分方程的数值解法实验室数学实验室所属课程名称微分方程数值解实验类型上机实验实验日期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\）的变化。
常微分方程实验报告
常微分方程实验报告《常微分方程》综合性实验实验报告实验班级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．常微分方程的数值解法除常系数线性微分方程可用特征根法求解、少数特殊方程可用初等积分法求解外，大部分微分方程无解析解，应用中主要依靠数值解法。
常微分方程的求解与定性分析实验报告
常微分方程的求解与定性分析实验报告Company Document number：WTUT-WT88Y-W8BBGB-BWYTT-19998常微分方程的求解与定性分析实验报告一、实验综述1、实验目的及要求●归纳和学习求解常微分方程(组)的基本原理和方法；●掌握解析、数值解法，并学会用图形观察解的形态和进行解的定性分析；●熟悉MATLAB软件关于微分方程求解的各种命令；●通过范例学习建立微分方程方面的数学模型以及求解全过程；●通过该实验的学习，使学生掌握微分方程(组)求解方法（解析法、欧拉法、梯度法、改进欧拉法等），对常微分方程的数值解法有一个初步了解，同时学会使用MATLAB软件求解微分方程的基本命令，学会建立微分方程方面的数学模型。
这对于学生深入理解微分、积分的数学概念，掌握数学的分析思维方法，熟悉处理大量的工程计算问题的方法是十分必要的。
2、实验仪器、设备或软件电脑、二、实验过程（实验步骤、记录、数据、分析）实验内容：根据实验内容和步骤，完成以下实验，要求写出实验报告（实验目的→问题→数学模型→算法与编程→计算结果→分析、检验和结论）1．求微分方程的解析解，并画出它们的图形。
y '= y + 2 x, y (0) = 1, 0< x <1；m=dsolve('Dy=y+2*x','y(0)=1','x')ezplot(m,[0 1])m =3*exp(x) - 2*x – 21．求微分方程⎪⎩⎪⎨⎧====-+]100[0)0(;0)0(01.03t uu u u u 的数值解，要求编写求解程序。
function dy=vdp1000(t,y)dy=zeros(2,1);dy(1)=y(2);dy(2)=-y(1)+*y(1)^3;[T,Y]=ode15s('vdp1000',[0 10],[0 0]);plot(T,Y(:,1),'-')3．Rossler 微分方程组：当固定参数b =2，c =4时，试讨论随参数a 由小到大变化（如 a ∈(0，)而方程解的变化情况，并且画出空间曲线图形，观察空间曲线是否形成混沌状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::[t,x]=ode45('rossler',t0,[0,0,0]);subplot(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:a=:a=:a=:a=:a=:结果分析：从图像可以看出，当a=0时，微分方程的解（x,y,z）收敛与（0，，）；当a=时，（x,y,z）仍收敛与（0，，），只是收敛速度减慢；当a=时，（x,y,z）已发散，周期性变化；随着a的增大，(x,y,z)接近其极限环的速度加快，空间曲线成混沌状。
大学数学实验四_常微分方程数值解
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 法和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,则超出运算的范围，发生溢出。
常微分方程上机实验报告
实验名称
常微分方程初值问题的数值解法
实验时间
2012年3月13日
姓名
XXX
班级
XXX
学号
XXX
成绩
一、实验目的及内容
1．MATLAB基本操作：矩阵相乘、矩阵LU分解及矩阵分解为另一个矩阵和其转置相乘、解矛盾方程、符号求解微分方程。
2．MATLAB数据调用及绘图。
二、相关背景知识介绍
（4）符号求解微分方程，求 。
2.MATLAB数据调用及绘图：
三．代码
1.(1)A=[3 4;5 6]
B=[0 8;1 6]
A*B
(2)A=[2 3;4 5]
[L U]=lu(A)
B=Chol(A)
(3)A=[2 8;3 4]
b=[6;2]
X=A\b
(4)dsolve('Dy=x','x')
2.
dsolve('Dy=x','y(0)=0','x')
1.MATLAB基本操作：矩阵Fra bibliotek乘、矩阵LU分解及矩阵分解为另一个矩阵和其转置相乘、解矛盾方程、符号求解微分方程。
（1）矩阵相乘：A=[3 4;5 6] B=[0 8;1 6]，求A*B。
（2）矩阵分解：A=[2 3;4 5]，求[L U]=lu(A)、B=Chol(A)。
（3）解方程：A=[2 8;3 4] b=[6;2]，求X=A\b。
ans =
-x-1+exp(x)*C1
>> dsolve('Dy=x','x')
ans =
1/2*x^2+C1
微分方程数值解上机报告-常微分方程初值问题数值解法
第一章 常微分方程初值问题数值解法实验一1、实验题目试用(a)欧拉格式(b)中点格式 (c)预报—校正格式 (d)经典四级四阶R-K 格式 编程计算方程：(]()⎪⎩⎪⎨⎧=∈+=112,1222y x y x x y dx dy 2、程序#include<iostream.h> #include<stdlib.h> #include<math.h> const int N=11;double fund(double x,double y); void Euler(double a,double h,double y[]); void Center(double a,double h,double y[]); void YuJiao(double a,double h,double y[]); void SjSj(double a,double h,double y[]); void Adams(double a,double h,double y[]); void main() { double a,b,h,y[N]; int i; char option;cout<<"请输入定义域区间[a,b]：\n";cin>>a>>b;cout<<"请输入初值y[0]：\n";cin>>y[0];h=(b-a)/(N-1);cout<<"a:欧拉法\nb:中心法\nc:预报-校正格式\n";cout<<"d:经典四级四阶R-K格式\n";cout<<"e:Adams预报-修正格式\nf:退出\n";label: cout<<"请选择算法：";cin>>option;switch(option){case 'a': Euler(a,h,y);break;case 'b': Center(a,h,y);break;case 'c': YuJiao(a,h,y);break;case 'd': YuJiao(a,h,y);break;case 'e': Adams(a,h,y);break;case 'f': exit(1);break;default : {cout<<"选择有错,重新选择！\n";goto label;}}cout<<"计算结果为：\n xn\t\tyn"<<endl;for(i=0;i<N;i++)cout<<" "<<a+i*h<<"\t\t"<<y[i]<<endl;}void Euler(double a,double h,double y[]){int i;for(i=1;i<N;i++){y[i]=y[i-1]+h*fund(a+(i-1)*h,y[i-1]);}}void Center(double a,double h,double y[]){int i;double w;for(i=1;i<N;i++){w=fund(a+(i-1)*h,y[i-1]);y[i]=y[i-1]+h*fund(a+(i-1)*h/2,w*h/2+y[i-1]);}}void YuJiao(double a,double h,double y[]){int i;double sun;double y_Begin[N]={y[0]};for(i=1;i<N;i++){y_Begin[i]=y_Begin[i-1]+h*fund(a+(i-1)*h,y_Begin[i-1]);}for(i=1;i<N;i++){while(fabs(y_Begin[i]-sun)>0.0001){sun=y[i-1]+h/2.0*(fund(a+(i-1)*h,y[i-1])+fund(a+i*h,y_Begin[i]));y_Begin[i]=sun;}y[i]=sun;}}void SjSj(double a,double h,double y[]){int i;double k1,k2,k3,k4;for(i=1;i<N;i++){k1=fund(a+(i-1)*h,y[i-1]);k2=fund(a+(i-1)*h+h/2,y[i-1]+h*k1/2);k3=fund(a+(i-1)*h+h/2,y[i-1]+h*k2/2);k4=fund(a+(i-1)*h+h,y[i-1]+h*k3);y[i]=y[i-1]+h/6*(k1+2*k2+2*k3+k4);}}void Adams(double a,double h,double y[]){int i;double me,x[N];double k1,k2,k3,k4;for(i=0;i<N;i++){x[i]=a+i*h;y[i]=y[0];}for(i=1;i<4;i++){k1=fund(a+(i-1)*h,y[i-1]);k2=fund(a+(i-1)*h+h/2,y[i-1]+h*k1/2);k3=fund(a+(i-1)*h+h/2,y[i-1]+h*k2/2);k4=fund(a+(i-1)*h+h,y[i-1]+h*k3);y[i]=y[i-1]+h/6*(k1+2*k2+2*k3+k4);}for(i=4;i<N;i++){me=y[i-1]+h/24*(55*fund(x[i-1],y[i-1])-59*fund(x[i-2],y[i-2])+37*fund(x[i-3],y[i-3])-9*fund(x[i-4],y[i-4]));while(fabs(me-y[i])>0.0001){y[i]=y[i-1]+h/24*(9*fund(x[i],me)+19*fund(x[i-1],y[i-1])-5*fund(x[i-2],y[i-2])+fund(x[i-3],y[i-3] ));me=y[i];}}}double fund(double x,double y){double s;s=y/(2*x)+x*x/2/y;return s;}3、试验结果及分析(i)运算结果(ii)结果分析在使用欧拉法数值求解过程中，其计算过程非常简单，即由y 可直接计算出1y ，由1y 可直接计算出,,2 y 无需用迭代方法求解任何方程，就可求出近似解ny 。
数学实验基础 实验报告(1)常微分方程
实验一 常微分方程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.4074实验一 常微分方程0.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.641500.20.40.60.81 1.2 1.4 1.6若x 上限增为1.58，1.60,则超出运算的范围，发生溢出。
常微分方程实验报告_闵安游_U201310085
⑵、3 步 4 阶 Adams 方法： ①、画图程序：
w = exp(i*linspace(0,2*pi)); z = 24*(w.^3 - w.^2)./(1 - 5*w + 19*w.^2 + 9*w.^3); plot(z)
②、绝对稳定性判断程序：
9 3 19 2 5 h (1 h ) (1 h) h 0 此时的特征方程是： 24 24 24 24
fprintf('误差为: %10.10f\n',m)
测试程序：
f = @(t,y)-2*y + 2*t + 1; F = @(t)exp(-2*t) + t; Threesteps_Adams1(0,2,f,F,0.1)
第 2 题： ⑴、2 步 3 阶 Adams 方法： ①、画图程序：
w = exp(i*linspace(0,2*pi)); z = 12*w.*(w - 1)./(5*w.*w + 8*w -1); plot(z)
t = a + i*h; y2 = y1 + h/2*(3*f1 - f0); plot(t,y2,'.'); f2 = f(t,y2); y1 = y2; f0 = f1; f1 = f2; end y = linspace(a,b); plot(y,F(y),'r'); hold off m = abs(y1-F(b)); fprintf('数值解: %10.10f\n',y1) fprintf('误差为:%10.10f\n',m)
' u2 24u1 51u2 9 cos t 1/ 3sin t
2013春数学实验基础 实验报告(1)常微分方程
实验一 常微分方程1. 分别用Euler 法和ode45解下列常微分方程并与解析解比较： (1)30,1)0(,<<=+='x y y x y编写Euler 法的matlab 函数，命名为euler.mfunction [t,y]=euler(odefun,tspan,y0,h) t=tspan(1):h:tspan(2); y(1)=y0;for i=1:length(t)-1y(i+1)=y(i)+h*feval(odefun,t(i),y(i)); end t=t'; y=y';下面比较三者的差别： % ode45odefun=inline('x+y','x','y'); [x1,y1]=ode45(odefun,[0,3],1); plot(x1,y1,'ko');pause hold on ; % Euler·¨[x2,y2]=euler(odefun,[0,3],1,0.05); plot(x2,y2,'r+');pause hold on ;% 解析解y0=dsolve('Dy=t+y','y(0)=1'); ezplot(y0,[0,3]);pause hold off ;legend('ode45','euler 法','解析解');实验一 常微分方程Euler 法只有一阶精度，所以实际应用效率比较差，而ode45的效果比较好，很接近真实值。
(2) 20.01()2sin(),(0)0,(0)1,05y y y t y y t ''''-+===<<先写M 文件ex1_2fun.mfunction f=ex1_2fun(t,y) f(1)=y(2);f(2)=0.01*y(2).^2-2*y(1)+sin(t); f=f(:); % ode45[t1,y1]=ode45(@ex1_2fun,[0,5],[0;1]); plot(t1,y1(:,1),'ko');% 解析解s=dsolve('D2y-0.01*(Dy)^2+2*y=sin(t)','y(0)=0','Dy(0)=1','t')s =[ empty sym ]%由此可知该微分方程无解析解2. 求一通过原点的曲线，它在(,)x y 处的切线斜率等于22,0 1.57.x y x +<<若x 上限增为1.58，1.60会发生什么？odefun=inline('2*x+y^2','x','y'); subplot(1,4,1);[x1,y1]=ode45(odefun,[0,1.57],0); plot(x1,y1,'r*'); title('上限1.57'); subplot(1,4,2);[x2,y2]=ode45(odefun,[0,1.58],0); plot(x2,y2,'bo'); title('上限1.58');实验一 常微分方程subplot(1,4,3);[x3,y3]=ode45(odefun,[0,1.6],0); plot(x3,y3,'k'); title('上限1.60'); subplot(1,4,4);plot(x1,y1,'r*');hold on ; plot(x2,y2,'bo');hold on ; plot(x3,y3,'k');hold off ;legend('上限1.57','上限1.58','上限1.60');上限1.5713上限1.5814上限1.6014结论：随着x 上界的增加，解趋于无穷大。
常微分方程实验报告.
常微分方程课程实验报告实验名称 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
常微分方程实验报告(一)实验名称 姓名： 学 号 班级：一、 实验目的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 =⎧⎪=⎨⎪=⎩。
微分方程数值解第四次上机报告
南京信息工程大学 实验（实习）报告实验课程 微分方程数值解 实验名称 第四次实验 实验日期 指导老师 王曰朋 专业 数学与应用数学 年级 姓名 学号 得分- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -实验目的：*编程五点三层格式求解两道二阶波动方程定解问题； *讨论不同网格划分得到的数值求解结果. 实验内容： 1. 编程内容对二阶波动方程定解问题（1）求解的程序见wave.m 文件，对定解问题（2）求解的程序见wave2.m 文件.2. 问题求解*先求解如下一道二阶波动方程定解问题（t 范围取有限）：⎪⎩⎪⎨⎧≤≤==≤≤==<<<<=--.10,0)0,(,)0,(,10,0),1(),0(10,10,u 2)3.0(400x x u e x u t t u t u t x u t x xx tt （1） 数值求解：用五点三层差分格式求解，先取空间步长02.0=∆x ,时间步长01.0=∆t ,故柯朗数15.0<=∆∆=xtc r ,满足稳定性条件.输出数值解三维图和末时刻解如图1.图1：疏网格)5.0(=r 时数值解三维图（左）和末时刻解图像（右）从图1可以看出数值解图像纹理粗糙，末时刻解还有微量波动，这时可以加密网格获得更好的解.为加密网格，再取空间步长002.0=∆x ,时间步长001.0=∆t .柯朗数不变，仍满足稳定性条件.输出数值解三维图和末时刻解如图2.从图2可以看出数值解三维图纹理相对图1更为平滑，且直观上看末时刻解比较稳定。
图2：较密网格)5.0(=r 时数值解三维图（左）和末时刻解图像（右） 再改动网格划分，取空间步长001.0=∆x ,时间步长001.0=∆t ,柯朗数 11≤=∆∆=xt c r ，仍满足稳定性条件. 输出数值解三维图和末时刻解如图3。
实验9 常微分方程
在图7.3中画出积分曲线，观察其吻合程度。
西安理工大学应用数学系
如若C=1时，画积分曲线的Matlab代码为 function jifenquxian clear; C=1; x=-4:0.01:4; for i=1:length(x) y(i)=-(exp(2*C - 2*x(i)) + 1)./(exp(2*C - 2*x(i)) - 1); plot(x(i),y(i),'r.'),pause(0.005) hold on end hold on 运行结果如图7.4.
西安理工大学应用数学系
（2）是二阶方程的初值问题，求特解，相应的Matlab代码为 clear; s=dsolve('D2y=sin(2*x)-y','y(0)=0','Dy(0)=1','x') y=simplify(s) %以最简形式显示s 运行结果为s =5/3*sin(x)-1/3*sin(2*x)， y = -1/3*(-5+2*cos(x))*sin(x) （3）是一阶方程组的初值问题，求特解，相应的Matlab代码为 clear; s=dsolve('Dx=x+y','Dy=y-x','x(0)=1','y(0)=1'); x=simplify(s.x) %s是一个结构 y=simplify(s.y) 运行结果为x =exp(t)*(cos(t) + sin(t)), y =exp(t)*(cos(t) - sin(t))
v f 1 f 2 ( x0 , y( x0 ))
常微分方程应用实验内容
实验一 MatLab 入门 实验目的： 实验内容与步骤：1.用MatLab 的命令方式计算表达式2.459.3*)43sin(2*)3ln(8÷+π的值.2. 生成矩阵A ,通过修改矩阵A 的第1行第3列元素为2,把矩阵扩充为4行5列，其中第4行第3列元素为7，其余列为0,删除第2列元素的操作生成矩阵B , 按行的逆顺序取A 的1,2,4,5列,且第四行为自然数1到4，生成新的矩阵C ;计算B *C ，B .* C , B +C ,2*B , B 和C 的行列式，B \C , B /C ,并找出A 中大于3的元素，且将其替换为1。
A = 0 1 0 2 1 3 4 6 8 4 9 7 3 2 23.在同一画面绘制0≤x ≤2*pi 范围内的sin(2x)、sinx 2、sin 2x 的图形。
4.在同一平面中的两个窗口分别用polar 绘制心形线)cos 1(3θρ-=和马鞍面z=x*x-y*y 的图形，并以不同的角度观察马鞍面。
5.在区域-3<x<3,-3<y<3绘制二元函数z=0.1*sin(x 2+y 2)的图形。
6.用ezplot 绘制函数e xy -sin(x+y)=0在[-3,3]上的图形。
7.用ezplot 绘制函数]6,0[,)cos 1(2)sin (2π∈⎩⎨⎧-=-=t t y t t x 的图形。
8用MatLab 软件求极限xx x cos 112ln 2lim 0x ---→9用MatLab 软件求下列函数的导数： （1）已知,)1()3(254+-+=x x x y 求1'=x y（2）已知x y x cos e =，求）（4y10 使用MatLab 软件求下列积分： （1）x xxx d sin cos 3⎰(2) x e x x xd 1e 2⎰+）（ (3)dx x e x ⎰202cos π（4）dx e x ⎰-10211 求解微分方程1)10()/log('==y x y y xy上机实验结果11. (1) 在windows 窗口，点击开始－>附件－>计算器，打开计算器，(2) 按如下顺序点击各键：(注：下面说明中各键之间用，隔开)选中十进制，弧度，C , 3, ln, *, (, 2, x^y, 8, ), =, MS, C, 3, /, 4, *, pi, =, sin, *, (, 3, ., 5, 9, x^y, 0, ., 5 , ), /, 4, ., 2, =, M+ , MR记录结果2.459.3*)43sin(2*)3ln(8÷+π＝ 。
常微分实习
目录1.实习的目的和任务2.实习要求3.实习地点4.主要仪器设备5.实习内容5.1 用不同格式对同一个初值问题的数值求解及其分析5.1.1求精确解5.1.2用欧拉法求解5.1.3用改进欧拉法求解5.1.4用4级4阶龙格—库塔法求解5.1.5 问题讨论与分析5.2 一个算法不同不长求解同一个初值问题及其分析6.结束语参考文献常微分方程课程实习1. 实习的目的和任务目的：通过课程实习能够应用MATLAB软来计算微分方程（组）的数值解；了解常微分方程数值解。
任务：通过具体的问题，利用MATLAB软件来计算问题的结果，分析问题的结论。
2. 实习要求能够从案例的自然语言描述中，抽象出其中的数学模型；能够熟练应用所学的数值解计算方法；能够熟练使用MATLAB软件；对常微分方程数值解有所认识，包括对不同算法有所认识和对步长有所认识。
3. 实习地点数学实验室、学生宿舍4. 主要仪器设备计算机、Microsoft Windows XP Matlab 7.05. 实习内容5.2 用不同格式对同一初值问题的数值求解及其分析用欧拉方法，改进欧拉方法，4阶龙格—库塔方法分别求下面微分方程的初值dy/dx=cos(x).^2-(tan(x))*y y(0)=2 x ∈(-2π，2π) 1.1 .1求精确解首先可以求得其精确解为：y=sin(x)*cos(x)+2*cos(x) 5.2.1 程序代码：x=-2*pi:0.1:2*pi;y=sin(x).*cos(x)+2*cos(x) plot(x,y,'b*-'); Data=[x',y' ]-8-6-4-22468-2.5-2-1.5-1-0.500.511.522.55.2.2用欧拉法求解 程序如下：建立函数文件cwfa1.mfunction [x,y]=cwfa1(fun,x_span,y0,h) x=x_span(1):h:x_span(2); y(1)=y0;for n=1:length(x)-1y(n+1)=y(n)+h*feval(fun,x(n),y(n)); end x=x';y=y';在MATLAB 输入以下程序：clear allfun=inline(' cos(x).^2-(tan(x))*y '); [x,y]=cwfa1(fun,[-2*pi,2*pi],2,0.1); [x,y]plot(x,y,'r*-')-8-6-4-202468-4-3-2-1012345.2.3用改进欧拉法求解： 程序如下：建立函数文件cwfa2.mfunction [x,y]=cwfa 2(fun,x_span,y0,h) x=x_span(1):h:x_span(2); y(1)=y0;for n=1:length(x)-1 k1=feval(fun,x(n),y(n)); y(n+1)=y(n)+h*k1; k2=feval(fun,x(n+1),y(n+1));y(n+1)=y(n)+h*(k1+k2)/2; end x=x';y=y';在MATLAB 输入以下程序：clear allfun=inline('cos(x).^2-(tan(x))*y'); [x,y]=cwfa2(fun,[-2*pi,2*pi],2,0.1); [x,y]plot(x,y,'y+-')-8-6-4-202468-3-2-112345.2.4 用4阶龙格—库塔求解 程序如下：建立函数文件cwfa3.mfunction [x,y]=cwfa3(fun,x_span,y0,h) x=x_span(1):h:x_span(2); y(1)=y0;for n=1:length(x)-1 k1=feval(fun,x(n),y(n));k2=feval(fun,x(n)+h/2,y(n)+h/2*k1); k3=feval(fun,x(n)+h/2,y(n)+h/2*k2); k4=feval(fun,x(n+1),y(n)+h*k3);y(n+1)=y(n)+h*(k1+2*k2+2*k3+k4)/6; end x=x';y=y';在MATLAB 输入以下程序： clear all;fun=inline(' cos(x).^2-(tan(x))*y'); [x,y]=cwfa3(fun,[-2*pi,2*pi],2,0.1); [x,y]plot(x,y, 'g+-')-8-6-4-22468-2.5-2-1.5-1-0.500.511.522.55.2.5 问题讨论与分析 由以上数值分析结果绘制表格：x=-2*pi:0.1:2*pi;y=sin(x).*cos(x)+2*cos(x) plot(x,y,'b*-'); hold on Data=[x',y']fun=inline(' cos(x).^2-(tan(x))*y '); [x,y1]=cwfa1(fun,[-2*pi,2*pi],2,0.1); plot(x,y1,'r*-')fun=inline('cos(x).^2-(tan(x))*y'); [x,y2]=cwfa2(fun,[-2*pi,2*pi],2,0.1); plot(x,y2,'y+-')fun=inline(' cos(x).^2-(tan(x))*y');[x,y3]=cwfa3(fun,[-2*pi,2*pi],2,0.1);plot(x,y3, 'g+-')D=[x,y',y1,y1-y',y2,y2-y',y3,y3-y']%利用矩阵D=-6.2832 2.0000 2.0000 -0.0000 2.0000 -0.0000 2.0000 -0.0000-6.1832 2.0893 2.1000 0.0107 2.0890 -0.0004 2.0893 -0.0000 -6.0832 2.1548 2.1779 0.0231 2.1541 -0.0008 2.1548 -0.0000 -5.9832 2.1930 2.2298 0.0368 2.1918 -0.0012 2.1930 -0.0000 -5.8832 2.2008 2.2521 0.0513 2.1991 -0.0017 2.2008 -0.0000 -5.7832 2.1759 2.2417 0.0658 2.1737 -0.0022 2.1759 -0.0000 -5.6832 2.1167 2.1963 0.0796 2.1139 -0.0028 2.1167 -0.0000 -5.5832 2.0224 2.1142 0.0917 2.0191 -0.0033 2.0224 -0.0000 -5.4832 1.8932 1.9946 0.1014 1.8894 -0.0038 1.8932 -0.0000 -5.3832 1.7301 1.8377 0.1076 1.7259 -0.0043 1.7301 -0.0000 -5.2832 1.5353 1.6448 0.1096 1.5306 -0.0047 1.5352 -0.0000 -5.1832 1.3114 1.4178 0.1064 1.3065 -0.0050 1.3114 -0.0000 -5.0832 1.0624 1.1598 0.0974 1.0573 -0.0051 1.0624 -0.0000 -4.9832 0.7927 0.8746 0.0819 0.7877 -0.0051 0.7927 -0.0000 -4.8832 0.5074 0.5667 0.0593 0.5027 -0.0048 0.5074 -0.0000 -4.7832 0.2120 0.2410 0.0290 0.2077 -0.0044 0.2119 -0.0002 -4.6832 -0.0876 -0.0984 -0.0108 -0.0834 0.0042 -0.0829 0.0047 -4.5832 -0.3855 -0.4350 -0.0495 -0.3672 0.0183 -0.3647 0.0208 -4.4832 -0.6757 -0.7681 -0.0925 -0.6439 0.0317 -0.6390 0.0367 -4.3832 -0.9525 -1.0922 -0.1397 -0.9080 0.0445 -0.9003 0.0522 -4.2832 -1.2107 -1.4014 -0.1907 -1.1541 0.0565 -1.1435 0.0672 -4.1832 -1.4455 -1.6903 -0.2449 -1.3776 0.0679 -1.3640 0.0815 -4.0832 -1.6528 -1.9539 -0.3011 -1.5744 0.0784 -1.5578 0.0950 -3.9832 -1.8294 -2.1877 -0.3583 -1.7413 0.0881 -1.7219 0.1075 -3.8832 -1.9729 -2.3881 -0.4153 -1.8760 0.0969 -1.8539 0.1190 -3.7832 -2.0817 -2.5525 -0.4708 -1.9771 0.1047 -1.9524 0.1293 -3.6832 -2.1555 -2.6790 -0.5235 -2.0441 0.1114 -2.0172 0.1383-3.5832 -2.1945 -2.7667 -0.5722 -2.0775 0.1171 -2.0486 0.1459 -3.4832 -2.2001 -2.8158 -0.6157 -2.0785 0.1215 -2.0480 0.1521 -3.3832 -2.1742 -2.8271 -0.6529 -2.0494 0.1248 -2.0175 0.1567 -3.2832 -2.1197 -2.8025 -0.6828 -1.9928 0.1269 -1.9599 0.1598 -3.1832 -2.0398 -2.7445 -0.7046 -1.9121 0.1277 -1.8786 0.1613 -3.0832 -1.9383 -2.6561 -0.7177 -1.8110 0.1273 -1.7772 0.1611 -2.9832 -1.8192 -2.5409 -0.7217 -1.6935 0.1257 -1.6598 0.1594 -2.8832 -1.6865 -2.4028 -0.7162 -1.5638 0.1228 -1.5305 0.1560 -2.7832 -1.5444 -2.2458 -0.7014 -1.4258 0.1187 -1.3933 0.1511 -2.6832 -1.3967 -2.0740 -0.6773 -1.2833 0.1134 -1.2519 0.1447 -2.5832 -1.2468 -1.8912 -0.6444 -1.1399 0.1070 -1.1100 0.1369 -2.4832 -1.0980 -1.7011 -0.6032 -0.9985 0.0995 -0.9703 0.1277 -2.3832 -0.9526 -1.5070 -0.5544 -0.8615 0.0911 -0.8354 0.1172 -2.2832 -0.8126 -1.3115 -0.4989 -0.7309 0.0817 -0.7071 0.1055 -2.1832 -0.6793 -1.1169 -0.4376 -0.6078 0.0715 -0.5865 0.0928 -2.0832 -0.5532 -0.9249 -0.3717 -0.4926 0.0607 -0.4741 0.0791 -1.9832 -0.4344 -0.7364 -0.3020 -0.3852 0.0492 -0.3697 0.0647 -1.8832 -0.3222 -0.5520 -0.2298 -0.2849 0.0373 -0.2726 0.0496 -1.7832 -0.2155 -0.3717 -0.1561 -0.1904 0.0251 -0.1815 0.0340 -1.6832 -0.1129 -0.1949 -0.0820 -0.1002 0.0127 -0.0948 0.0181 -1.5832 -0.0124 -0.0210 -0.0086 -0.0142 -0.0018 -0.0111 0.0013 -1.4832 0.0878 0.1482 0.0603 0.1003 0.0125 0.0787 -0.0092 -1.3832 0.1898 0.3176 0.1279 0.2163 0.0265 0.1702 -0.0196 -1.2832 0.2953 0.4884 0.1931 0.3354 0.0401 0.2655 -0.0298 -1.1832 0.4060 0.6616 0.2556 0.4593 0.0532 0.3663 -0.0397 -1.0832 0.5231 0.8380 0.3148 0.5889 0.0658 0.4739 -0.0492 -0.9832 0.6474 1.0179 0.3706 0.7249 0.0776 0.5891 -0.0582 -0.8832 0.7789 1.2015 0.4226 0.8675 0.0886 0.7123 -0.0667 -0.7832 0.9173 1.3880 0.4707 1.0160 0.0986 0.8429 -0.0744 -0.6832 1.0615 1.5764 0.5149 1.1692 0.1077 0.9801 -0.0815 -0.5832 1.2098 1.7649 0.5552 1.3254 0.1157 1.1221 -0.0877-0.4832 1.3596 1.9510 0.5914 1.4821 0.1225 1.2666 -0.0930 -0.3832 1.5082 2.1318 0.6236 1.6362 0.1280 1.4108 -0.0974 -0.2832 1.6521 2.3038 0.6517 1.7843 0.1323 1.5512 -0.1008 -0.1832 1.7874 2.4630 0.6756 1.9226 0.1352 1.6842 -0.1033 -0.0832 1.9103 2.6053 0.6950 2.0469 0.1367 1.8056 -0.1047 0.0168 2.0165 2.7263 0.7098 2.1533 0.1368 1.9115 -0.1050 0.1168 2.1021 2.8217 0.7196 2.2376 0.1355 1.9978 -0.1043 0.2168 2.1633 2.8873 0.7240 2.2960 0.1328 2.0607 -0.1026 0.3168 2.1965 2.9190 0.7225 2.3252 0.1287 2.0967 -0.0998 0.4168 2.1990 2.9136 0.7147 2.3223 0.1233 2.1029 -0.0960 0.5168 2.1684 2.8682 0.6998 2.2850 0.1166 2.0771 -0.0913 0.6168 2.1033 2.7808 0.6775 2.2120 0.1087 2.0176 -0.0857 0.7168 2.0031 2.6501 0.6470 2.1028 0.0997 1.9239 -0.0792 0.8168 1.8681 2.4760 0.6079 1.9578 0.0897 1.7962 -0.07190.9168 1.6995 2.2592 0.5597 1.7784 0.0788 1.6356 -0.06391.0168 1.49962.0015 0.5019 1.5667 0.0672 1.4443 -0.0553 1.1168 1.2712 1.7056 0.4344 1.3261 0.0549 1.2252 -0.0461 1.2168 1.0184 1.3753 0.3569 1.0606 0.0422 0.9820 -0.0364 1.3168 0.7457 1.0152 0.2695 0.7749 0.0292 0.7193 -0.0264 1.4168 0.4583 0.6304 0.1721 0.4745 0.0162 0.4422 -0.0162 1.5168 0.1618 0.2266 0.0648 0.1646 0.0028 0.1559 -0.0059 1.6168 -0.1380 -0.1925 -0.0545 -0.1392 -0.0012 -0.0725 0.0655 1.7168 -0.4350 -0.6102 -0.1753 -0.4391 -0.0042 -0.2278 0.2071 1.8168 -0.7233 -1.0230 -0.2997 -0.7309 -0.0076 -0.3766 0.34671.9168 -0.9974 -1.4245 -0.4271 -1.0086 -0.0112 -0.5146 0.48272.0168 -1.2519 -1.8081 -0.5562 -1.2669 -0.0150 -0.6379 0.6140 2.1168 -1.4824 -2.1677 -0.6853 -1.5012 -0.0188 -0.7432 0.7391 2.2168 -1.6847 -2.4974 -0.8127 -1.7072 -0.0225 -0.8278 0.8569 2.3168 -1.8559 -2.7924 -0.9366 -1.8819 -0.0261 -0.8898 0.9661 2.4168 -1.9936 -3.0485 -1.0549 -2.0230 -0.0294 -0.9280 1.0656 2.5168 -2.0966 -3.2624 -1.1658 -2.1291 -0.0324 -0.9421 1.15452.6168 -2.1645 -3.4319 -1.2675 -2.1996 -0.0352 -0.9327 1.2318 2.7168 -2.1978 -3.5557 -1.3580 -2.2353 -0.0375 -0.9009 1.2969 2.8168 -2.1979 -3.6335 -1.4357 -2.2374 -0.0395 -0.8489 1.34892.9168 -2.1670 -3.6661 -1.4991 -2.2080 -0.0411 -0.7794 1.38753.0168 -2.1079 -3.6548 -1.5469 -2.1501 -0.0422 -0.6957 1.4123 3.1168 -2.0242 -3.6022 -1.5781 -2.0670 -0.0428 -0.6012 1.4229 3.2168 -1.9194 -3.5112 -1.5918 -1.9624 -0.0430 -0.5001 1.4193 3.3168 -1.7977 -3.3853 -1.5876 -1.8404 -0.0427 -0.3962 1.4016 3.4168 -1.6632 -3.2284 -1.5652 -1.7052 -0.0420 -0.2934 1.3698 3.5168 -1.5199 -3.0446 -1.5248 -1.5607 -0.0409 -0.1956 1.3243 3.6168 -1.3715 -2.8382 -1.4666 -1.4108 -0.0393 -0.1059 1.2656 3.7168 -1.2217 -2.6131 -1.3914 -1.2590 -0.0373 -0.0274 1.1943 3.8168 -1.0732 -2.3732 -1.3000 -1.1082 -0.0350 0.0378 1.11103.9168 -0.9286 -2.1223 -1.1936 -0.9609 -0.0322 0.0880 1.01664.0168 -0.7897 -1.8633 -1.0736 -0.8189 -0.0292 0.1224 0.9121 4.1168 -0.6576 -1.5990 -0.9414 -0.6834 -0.0259 0.1409 0.7985 4.2168 -0.5327 -1.3316 -0.7988 -0.5550 -0.0222 0.1441 0.6769 4.3168 -0.4151 -1.0626 -0.6475 -0.4335 -0.0184 0.1334 0.5485 4.4168 -0.3039 -0.7933 -0.4894 -0.3182 -0.0143 0.1107 0.4146 4.5168 -0.1980 -0.5243 -0.3263 -0.2080 -0.0100 0.0785 0.2766 4.6168 -0.0959 -0.2559 -0.1600 -0.1014 -0.0056 0.0399 0.1358 4.7168 0.0044 0.0119 0.0075 0.0116 0.0072 0.0030 -0.0014 4.8168 0.1048 0.2818 0.1770 0.2744 0.1696 0.0723 -0.03254.9168 0.2072 0.5518 0.3445 0.5375 0.3303 0.1439 -0.06335.0168 0.3135 0.8220 0.5085 0.8011 0.4876 0.2200 -0.0935 5.1168 0.4252 1.0927 0.6674 1.0652 0.6400 0.3025 -0.1228 5.2168 0.5435 1.3634 0.8199 1.3295 0.7860 0.3927 -0.1508 5.3168 0.6690 1.6337 0.9648 1.5931 0.9242 0.4916 -0.1773 5.4168 0.8017 1.9026 1.1009 1.8548 1.0531 0.5997 -0.2021 5.5168 0.9412 2.1684 1.2272 2.1126 1.1714 0.7164 -0.2248 5.6168 1.0862 2.4290 1.3428 2.3643 1.2780 0.8410 -0.24535.7168 1.2349 2.6818 1.4469 2.6068 1.3719 0.9716 -0.2633 5.8168 1.3848 2.9236 1.5388 2.8367 1.4519 1.1061 -0.2787 5.9168 1.5328 3.1505 1.6177 3.0503 1.5175 1.2415 -0.29136.0168 1.6755 3.3586 1.6831 3.2433 1.5678 1.3745 -0.3010 6.1168 1.8091 3.5433 1.7342 3.4115 1.6024 1.5013 -0.30776.2168 1.9294 3.7000 1.7706 3.5504 1.6210 1.6181 -0.3113-8-6-4-202468-4-3-2-101234在matlab 中输入以下代码：subplot(2,2,1); plot(x,y,'b*-'); hold onplot(x,y1,'r*-') subplot(2,2,2); plot(x,y,'b*-'); plot(x,y2,'r*-') subplot(2,2,3); plot(x,y,'b*-'); plot(x,y3,'r*-') subplot(2,2,4); plot(x,y,'b*-'); plot(x,y,'b*-') plot(x,y1,'r*-')plot(x,y2,'y*-') plot(x,y3,'g*-')-10-50510-4-2024-10-50510-4-2024-10-50510-4-2024-10-50510-4-2024第一幅是精确解与欧拉方法的对比，第二幅是精确解与改进欧拉方法的对比，第三幅则是精确解与龙格-库塔法德比较用excel 分别求出他们对精确解的最小二乘估计：CORREL(y,y1)=0.976423,CORREL(y,y2)=0.975323和CORREL(y,y3)=0..95654(CORREL 是excel 的最小二乘估计函数函数)由最小二乘结果可以看出欧拉法误差最大，而改进欧拉和龙格—库塔方法误差相对较小，但是从图像的拟合程度上来说，改进的欧拉方法是最理想的。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
（1）矩阵相乘：A=[3 4;5 6] B=[0 8;1 6]，求A*B。
（2）矩阵分解：A=[2 3;4 5]，求[L U]=lu(A)、B=Chol(A)。
（3）解方程：A=[2 8;3 4] b=[6;2]，求X=A\b。
（4）符号求解微分方程，求 。
2.MATLAB数据调用及绘图：
三．代码
1.(1)A=[3 4;5 6]
B=[0 8;1 6]
A*B
(2)A=[2 3;4 5]
[L U]=lu(A)
B=Chol(A)
(3)A=[2 8;3 4]
b=[6;2]
X=A\b
(4)dsolve('Dy=x','x')
2.
dsolve('Dy=x','y(0)=0','x')
ans =
-x-1+exp(x)*C1
>> dsolve('Dy=x','x')
ans =
1/2*x^2+C1
2.
五．计算结果的分析
在三种对常微分方程的数值解法单步法的求解过程中，能清晰的比较出三种方法的精度大小，Euler法的代数精确度为一阶，改进的Euler法具有二阶精度，而经典的R-K法具有四阶精度，所以在Matlab的图像中能清楚地表现出来。
教师评语
指导教师：年月日
常微分方程及其数值解实验报告
实验名称
常微分方程初值问题的数值解法
实验时间
2012年3月13日
姓名
XXX
班级
XXX
学号
XXX
成绩
一、实验目的及内容
1．MATLAB基本操作：矩阵相乘、矩阵LU分解及矩阵分解为另一个矩阵和其转置相乘、解矛盾方程、符号求解微分方程。
2．MATLAB数据调用及绘图。
二、相关背景知识介绍U =2.0000 3.0000
0 3.5000
A=[2 3;4 5]
B=Chol(A)
A =
2 3
1 5
B =
1.4142 2.1213
0 0.7071
3.
A=[2 8;3 4]
B=[6;2]
X=A\b
A =
2 8
3 4
b =
6
2
x =
-0.5000
0.8750
4.
dsolve('Dy=y+x','x')
x=[0:0.03:2]
y=1/2*x.^2
plot(x,y,'r')
四.数值结果
1.
A=[3 4;5，6]
B=[0 8;1 6]
A*B
A =
2 3
4 5
B =
0 8
1 6
ans =
3 34
5 62
2.
A=[2 3;4 5]
[L U]=lu(A)
A =
2 3
1 5
L =
1.0000 0
0.5000 1.0000