当前位置：文档之家 › 常微分方程实验报告

常微分方程实验报告

  • 格式：doc
  • 大小：201.00 KB
  • 文档页数：6

下载文档原格式

  / 6

合集下载

常微分方程 课程实习报告

常微分方程 课程实习报告

福建农林大学计算机与信息学院（数学类课程）课程实习报告课程名称：常微分方程课程实习实习题目：常微分方程数值求解问题的实习姓名：上官火玉系：信息与计算科学专业：信息与计算科学年级：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由上表和图可以看出欧拉法误差最大，而改进欧拉和龙格—库塔方法误差相对较小，并且龙格—库塔方法误差最小且大部分值都跟精确值相同。

常微分方程数值解实验报告

常微分方程数值解实验报告

常微分方程数值解实验报告学院：数学与信息科学专业：信息与计算科学姓名：郑思义学号：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 法的精度更高。

常微分方程实验报告

常微分方程实验报告

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

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

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

对于一阶常微分方程，形如＼（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．常微分方程的数值解法除常系数线性微分方程可用特征根法求解、少数特殊方程可用初等积分法求解外，大部分微分方程无解析解，应用中主要依靠数值解法。

大学数学实验四_常微分方程数值解

大学数学实验四_常微分方程数值解
【实验内容】
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。

先随便猜测k的一个试验值,程序中令k=12.由Numerov算法根据本题的条件，kn+1=kn=kn-1=k,s=0,得到yn+2,yn+1，yn间的迭代公式令con=（k＊h）^2/12yn+2=2*(1-5＊con）＊yn+1/（1+con）—yn3自己给定φ的初始条件，然后利用公式得到边界值φ(1)4。

然后以小的步长dk增加k值，这里令dk=1,每当φ（1)改变符号时,就将步长减半后倒退回来重复5.当步长小于所要求的容许误差时终止程序，此时的k值即为所求。

（二）弦割法1.随便猜测两个k值,这里令k0=1,k1=22.自己给定φ的初始条件，对两个k值分别利用上述公式进行迭代，得到边界值y1(1）和y2（1）。

3.比较y1（1)和y2(1)的绝对值大小。

若绝对值大,说明对应的k值距离本征值距离较远。

4.将（k0+k1）/2赋给k2，边界值绝对值小的对应的k值保持不变，边界值绝对值大的对应k值重新定位k2的值.5.重复进行实验，当y1(1)和y(2）的差的绝对值小于容许误差时终止程序.此时k1的值即为所求。

当搜索法和弦割法大致求出了一个本征值后，利用打靶法，调整k值再度进行搜索，得到多个本征值，绘出其中一个本征值对应的函数图像，观察其性质。

三、程序实现1.搜索法subroutine add(t,y0，y1） !利用子程序表示函数值的迭代implicit nonereal(8）：：t,h，con，y0,y1，y2integer:：i，nn=10000h=1。

0/ncon=(t＊h)＊*2/12do i=1,n-1y2=2＊（1—5*con）＊y1/（1+con）—y0 ！利用Numerov算法,得到迭代公式y0=y1 !向前迭代y1=y2end doreturnend subroutine addprogram zy3implicit nonereal(8)：：diffk，dk，yold,k，b0,b1integer::sb0=0.01 !取初始值，根据题目条件，令y0=y1,来保证x=0的位置导数为0 b1=0。

数学实验基础 实验报告(1)常微分方程

数学实验基础 实验报告(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

常微分方程实验报告_闵安游_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

计算方法实验报告-常微分方程的数值解法

计算方法实验报告-常微分方程的数值解法
15.0000 50.1112 50.1112
15.1000 50.1224 50.1224
15.2000 50.1332 50.1332
15.3000 50.1436 50.1436
15.4000 50.1535 50.1535
disp([T',V',df']);
%画图观察效果
figure;
plot(T,df,'k*',T,V,'--r');
grid on;
title('速度v与时间t的函数曲线');
legend('准确值','四阶经典R-K法');
结果如下所示：
步长四阶经典R-K法准确值
0 0 0
0.1000 0.9799 0.9799
5.9000 41.1708 41.1708
6.0000 41.4919 41.4919
6.1000 41.8028 41.8028
6.2000 42.1039 42.1039
6.3000 42.3954 42.3954
6.4000 42.6775 42.6775
6.5000 42.9504 42.9504
11.5000 49.2667 49.2667
11.6000 49.3097 49.3097
11.7000 49.3511 49.3511
11.8000 49.3909 49.3909
11.9000 49.4292 49.4292
12.0000 49.4661 49.4661
12.1000 49.5016 49.5016
（2）某跳伞者在t=0时刻从飞机上跳出，假设初始时刻的垂直速度为0，且跳伞者垂直下落。已知空气阻力为F=cv2，其中c为常数，v为垂直速度，向下方方向为正。写出此跳伞者的速度满足的微分方程；若此跳伞者的质量为M=70kg，且已知c=0.27kg/m，利用四阶龙格-库塔公式计算t<=20s的速度（取h=0.1s）

数值分析实验报告之常微分方程数值解

数值分析实验报告之常微分方程数值解

数学与计算科学学院实验报告
实验项目名称常微分方程数值解
所属课程名称数值方法B
实验类型验证
实验日期2013.11.11
班级
学号
姓名
成绩
图1 h=0.1时三个方法走势图
图2 h=0.05时三个方法走势图
图4 h=0.05时三个方法走势图
附录1：源程序
附录2：实验报告填写说明
1．实验项目名称：要求与实验教学大纲一致。

2．实验目的：目的要明确，要抓住重点，符合实验教学大纲要求。

3．实验原理：简要说明本实验项目所涉及的理论知识。

4．实验环境：实验用的软、硬件环境。

5．实验方案（思路、步骤和方法等）：这是实验报告极其重要的内容。

概括整个实验过程。

对于验证性实验，要写明依据何种原理、操作方法进行实验，要写明需要经过哪几个步骤来实现其操作。

对于设计性和综合性实验，在上述内容基础上还应该画出流程图、设
计思路和设计方法，再配以相应的文字说明。

对于创新性实验，还应注明其创新点、特色。

6．实验过程（实验中涉及的记录、数据、分析）：写明具体实验方案的具体实施步骤，包括实验过程中的记录、数据和相应的分析。

7．实验结论（结果）：根据实验过程中得到的结果，做出结论。

8．实验小结：本次实验心得体会、思考和建议。

9．指导教师评语及成绩：指导教师依据学生的实际报告内容，给出本次实验报告的评价。

常微分方程实验报告.

常微分方程实验报告.

常微分方程课程实验报告实验名称 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

常微分方程实验报告(一)实验名称 姓名： 学 号 班级：一、 实验目的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)2. 实习要求 (1)3. 实习地点 (1)4. 主要仪器设备 (1)5. 实习内容 (1)5.1 用不同格式对同一个初值问题的数值求解及其分析 (1)5.1.1求精确解 (1)5.1.2用欧拉法求解 (6)5.1.3用改进欧拉法求解 (9)5.1.4用4级4阶龙格—库塔法求解 (12)5.1.5 问题讨论与分析 (15)5.2 一个算法不同不长求解同一个初值问题及其分析 (18)6. 结束语 (29)参考文献 (29)常微分方程课程实习1. 实习的目的和任务目的：通过课程实习能够应用MATLAB 软来计算微分方程（组）的数值解；了解常微分方程数值解。

任务：通过具体的问题，利用MATLAB 软件来计算问题的结果，分析问题的结论。

2. 实习要求能够从案例的自然语言描述中，抽象出其中的数学模型；能够熟练应用所学的数值解计算方法；能够熟练使用MATLAB 软件；对常微分方程数值解有所认识，包括对不同算法有所认识和对步长有所认识。

3. 实习地点数学实验室、学生宿舍 4. 主要仪器设备计算机、Microsoft Windows 7 Matlab 7.0 5. 实习内容5.1 用欧拉方法，改进欧拉方法，4阶龙格—库塔方法分别求下面微分方程的初值dy/dx=y*cos(x+2) y(-2)=1 x ∈[-2，0] 5.1.1求精确解 ①变量分离方程情形:形如()*()dyf xg y dx=的方程,这里(),()f x g y 分别是,x y 的连续函数.如果()0g y ≠,我们可将方程改写成()()dyf x dxg y =,这样,变量就”分离”开来了,两边同时积分即可:(),()dyf x dx c cg y =+⎰⎰为任意常数. ②常数变易法:一阶线性微分方程()*()dyf x yg x dx =+,其中(),()f x g x 在考虑区 间上是的连续函数.可先解出方程()*dyf x y dx=的解,这是属于变量分离方程情形,可解得:*exp(())y c f x dx =⎰,这里c 是任意常数.然后将变c易为x 的待定函数()c x ,令()*exp(())y c x f x dx =⎰,将其代入原方程可得:()*exp(())()*()*exp(())()*()*exp(())()dy dc x f x dx c x f x f x dx f x c x f x dx g x dx dx=+=+⎰⎰⎰所以可解得()()*exp(())*1c x f x f x dx dx c =-⎰⎰,这里1c 是任意常数.将()()*exp(())1c x f x f x dx dx c =-+⎰⎰代入()*exp(())y c x f x dx =⎰可得原方程的通解:(()*exp(())1)*exp(()),1y f x f x dx dx c f x dx c =-+⎰⎰⎰为任意常数.③恰当微分方程情形:形如(,)(,)0m x y dx n x y dy +=的一阶微分方程,这里 假设(,),(,)m x y n x y 在某矩形域内是的连续函数,且具有连续的一阶偏导数. 若m ny x∂∂=∂∂,则为恰当微分方程.判断为恰当微分方程后,则可用如下解法： 设(,)u x y c =是原方程的解，则um x∂=∂，所以设(,)()u m x y dx v y =+⎰， 则((,)())((,))()m x y dx v y m x y dx udv y n y y y dy ∂+∂∂==+∂∂∂⎰⎰＝，所以((,))()m x y dx dv y n ydy ∂-=∂⎰，由此()[(,)]v y n m x y dx dy y∂=-∂⎰⎰,由此可解得(,)[(,)]u m x y dx n m x y dx dy y ∂=+-∂⎰⎰⎰，所以原方程的通解为(,)[(,)],m x y dx n m x y dx dy c c y∂+-=∂⎰⎰⎰为任意常数。

  1. 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

《常微分方程》综合性实验实验报告实验班级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．常微分方程的数值解法除常系数线性微分方程可用特征根法求解、少数特殊方程可用初等积分法求解外，大部分微分方程无解析解，应用中主要依靠数值解法。

考虑一阶常微分方程初值问题000()(,()),()f y t f t y t t t t y t y '=<<⎧⎪⎨=⎪⎩ 其中12(,,,)T m y y y y = ，12(,,,)T m f f f f = ，010200(,,,)T m y y y y = 。

所谓数值解法，就是寻求()y t 在一系列离散节点01n f t t t t <<<≤ 上的近似值,0,1,,k y k n = ，称1k k k h t t +=-为步长，通常取为常量h 。

最简单的数值解法是Euler 法。

Euler 法的思路很简单：在节点处用差商近似代替导数1()()()k k k y t y t y t h+-'≈这样导出计算公式 1(,),0,1,2,k k k k y y hf t y k +=+=它能求解各种形式的微分方程。

Euler 法也称折线法。

Euler 法只有一阶精度，改进方法有二阶Runge-Kutta 法、四阶Runge-kutta 法、五阶Runge-kutta-Felhberg 法和线性多步法等，这些方法可用于解高阶常微分方程（组）初值问题。

边值问题采用不同方法，如差分法、有限元法等。

数值解法的主要据点是它缺乏物理意义。

3．用MA TLAB 数值求解微分方程的命令MATLAB 中主要用函数ode45，ode23，ode15s 求常微分方程的数值解。

其中ode45是最常用的求解微分方程数值解的命令，对于刚性方程组不宜采用。

ode23与ode45类似，只是精度低一些。

ode15s 用来求解刚性方程组。

这三个函数的调用格式基本相同，下面介绍ode45的调用格式。

[tout,yout]=ode45(‘yprime ’,[t0,tf],y0)此命令是采用变步长四阶Runge-kutta 法和五阶Runge-kutta-Felhberg 法求数值解，yprime 是表示f(t,y)的M 文件名，t0表示自变量的初始值，tf 表示自变量的终值，y0表示初始向量值。

输出向量tout 表示节点0,1(,,)T n t t t ，输出矩阵yout 表示数值解，每一列对应y 的一个分量。

若无输出参数，则自动作出图形。

四、实验方法与步骤1求下列常微分方程的通解(1)dx at dt =- 22(2)(s i n c o s )d y d y a b x x d x d x +=+ 2(3)(1)dy y y dx =- （4）22sin d y y dx= 求通解相应的MATLAB 代码为>> x1=dsolve('Dx=-a*t')>> y2=dsolve('D2y+a*Dy=b*(sin(x)+cos(x))','x')>> y3=dsolve('Dy=y^2*(1-y)','x')>> y4=dsolve('D2y=sin(y)','x')运行后屏幕显示常微分方程（1）、（2）、（3）、（4）的通解x1、y1、y2、y3、y4依次为： x1 =-1/2*a*t^2+C1y2 =b*(-a*cos(x)+a*sin(x)-sin(x)-cos(x))/(a^2+1)+C1+C2*exp(-a*x)y3= x+1/y-log(y)+log(-1+y)+C1=0y4=[ Int(1/(-2*cos(a)+C1)^(1/2),a=``..y)-x-C2=0,-Int(1/(-2*cos(a)+C1)^(1/2),a=``..y)-x-C2=0]2求下列常微分方程在给定初始条件下的特解。

22(1)1,(0)0dy y y dx ⎛⎫+== ⎪⎝⎭2220()()()(2),(0)1,x ad f x d f x d f x a f dx dx dx π===-=求特解相应的MATLAB 代码为>> y=dsolve('(Dy)^2+y^2=1','y(0)=0','x'),>> f=dsolve('D2f=-a^2*(Df)','f(0)=1,Df(pi/a)=0','x'),运行后屏幕显示常微分方程（1）、（2）在给定初始条件下的特解依次如下：y =[ sin(x)][ -sin(x)]f =13求解微分方程1,(0)1dy y t y dt=-++=的解析解和数值解，并进行比较。

相应的MA TLAB 代码为>> y=dsolve('Dy=-y+t+1','y(0)=1'),运行结果即方程的解析解为。

y =t+exp(-t)下面再求其数值解，先编写M 文件fun1.m 。

function f=fun1(t,y)f=-y+t+1；再运行相应的MATLAB 代码>> clear;close;t=0:0.1:1;>> y=t+exp(-t);plot(t,y); %画解析解的图形>> hold on; %保留已经画好的图形，如果下面再画图，两个图形合并在一起 >> [t,y]=ode45('fun1',[0,1],1);>> plot(t,y,'ro'); %画数值解图形，用红色小圈画>> xlabel('t'),ylabel('y');运行结果画图分析解析解和数值解吻合情况4.综合实验（数学建模）：求解如下问题之一：1）．小船从河边点o 处出发驶向以对岸（两岸为平行直线），设船速为a ，船行方向始终与河岸垂直，又设河宽为h ，河中任一点处的水流速度与该点到两岸距离的乘积成正比（比例系数为k ），求小船的航行路线及船到达对岸的位置。

2）．（目标的跟踪问题）设位于坐标原点的甲舰向位于x 轴上点)0,1(A 处的乙舰发射制导导弹，导弹头始终对准乙舰．如果乙舰以最大的速度0v （0v 是常数）沿平行于y 轴的直线行驶，导弹的速度是50v ，求导弹运行的曲线方程及乙舰行驶多远时，它将被导弹击中？ （说明：以上两题选做其中之一，或者自选其它题目）一． 问题：小船从河边点o 处出发驶向以对岸（两岸为平行直线），设船速为a ，船行方向始终与河岸垂直，又设河宽为h ，河中任一点处的水流速度与该点到两岸距离的乘积成正比（比例系数为k ），求小船的航行路线及船到达对岸的位置。

二． 实验目的：1．常微分方程的解了解析解和数值解。

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

三． 实验内容及要求：以船的起点为原点，水流方向为y 轴正向，船的航行方向为x 方向建立坐标系。

设在t 时刻船在x 轴方向上运行了x ，在y 轴方向上运行了y ，在很小的时间dt 内，可以认为船在y 轴方向作均速v 运动，移动了dy=vdt;在x 轴方向运动了dx=adt;其中v=k*x*(h-x),由此可建立微分方程组：**()(0)dy k x h x h dt t dx a a dt⎧=-⎪⎪≤≤⎨⎪=⎪⎩ 当t=0时，x 和y 都为0，将方程组中的t 消去得下微分方程：**()dy k x h x dx a=- ()0t h ≤≤ ()00y =输入相应的MA TLAB 代码：>> y=dsolve('Dy=(k/a)*x*(h-x)','y(0)=0','x')运行结果即方程的解析解为y =-1/3*k/a*x^3+1/2*k/a*h*x^2这也是小船的运动轨迹方程。

当x=h 时，船已运行到对岸，此时y 的值为36kh a ,则船在坐标中的位置是（h, 36kh a）. 下面就110.002**k s m --=, 2/a m s =,h=50m.画图看一下小船的运行曲线。

此时方程化为：3211**300040y x x =-+ (050x ≤≤) 下面画解析式的图其数值解图，先编写M 文件function f=fun1(x,y)f=(-1/3000)*x^3+1/40*x^2;function f=fun2(x,y)f=1/1000*x*(50-x);再运动相应的MATLAB 代码>> clear;close;x=0:1:50;fplot('fun1',[0,50]); %画解析解的图形hold on; %保留已经画好的图形，如果下面再画图，两个图形合并在一起[x,y]=ode45('fun2',[0,50],0);plot(x,y,'g*'); %画数值解图形，用绿色星号xlabel('x'),ylabel('y');运行结果画图四．结果分析：两曲线比较好地吻合在一起，此也为小船的运行曲线，由图可知小船的坐标为（50,125/6），也即是小船运动到对岸时已顺不流方向运动了20.8米。