MATLAB在求解常微分方程中的应用.docx
《MATLAB在求解常微分方程中的应用.docx》由会员分享,可在线阅读,更多相关《MATLAB在求解常微分方程中的应用.docx(19页珍藏版)》请在冰豆网上搜索。
MATLAB在求解常微分方程中的应用
《MATLA语言》课程论文
MATLAE在求解常微分方程中的应用
姓名:
袁学婷
学号:
12010245278
专业:
通信工程
班级:
通信一班
指导老师:
汤全武
学院:
物理电气信息学院
完成日期:
2012年12月10日
MATLAB在求解常微分方程中的应用
(袁学婷120102452782010通信班)
【摘要】MATLAB成为许多学科的解题工具,将MATLAB融入其它课程的学习中,可以大大提高运
算效率和准确性。
随着计算机的普及和国民整体素质的提高,科学计算将会更加的普及。
MATLAB在矩
阵及数值计算、多项式和线形代数、符号数学的基本方法等方面都有较好的应用。
本文概括地介绍了
MATLAB在求解常微分方程中的应用。
【关键词】MATLAB求解常微分方程解析解数值解
一、问题的提出
自20世纪80年代以来,出现了多种科学计算语言,亦称数学软件,比较流行的有MATLABMathematica、Maple等。
因为他们具有功能强、效率高、简单易学等特点,在在许多领域等到广泛应用。
MATLA便是一种影响大、流行广的科学计算语言。
MATLA的语法规则简单,更加贴近人的思维方式。
MATLA是英文MATrixLABoratory(矩阵实验室)的缩写。
自1984年由美国MathWorks公司推向市场以来,得到了广泛的应用和发展。
在欧美各高等院校MATLA已经成为线性代数、自动控制理论、数字信号处理、时间序列分析、动态系统仿真、图像处理等诸多课程的基本教学工具,成为大学生、硕士生以及博士生必须掌握的基本技能。
在设计研究单位和工业部门,MATLA已被广泛的应用于研究和解决各种具体的工程问题。
近年来,MATLAB
在我国也开始流行,应用MATLAB勺单位和个人急剧增加。
可以预见,MATLAB将在我国科
学研究和工程应用中发挥越来越大的作用。
现在简单的介绍一下MATLABE解微分方程中的
应用。
众所周知,只有对一些典型的常微分方程,才能求出它们的一般解表达式并用初始条件确定表达式中的任意常数。
所以能够求解的微分方程是十分有限的,特别是高阶方程
(组)•这就要求我们必须研究微分方程(组)的解法,既要研究微分方程(组)的解析解法(精确解),更要研究微分方程(组)的数值解法(近似解)•
对微分方程(组)的解析解法(精确解),Matlab有专门的函数可以用,本文将作一定的介绍.
二、常微分方程的概述
1、微分方程的概念
未知的函数以及它的某些阶的导数连同自变量都由一已知方程联系在一起的方程称为微分方程。
如果未知函数是一元函数,称为常微分方程。
常微分方程的一般形式为
F(t,y,y',y",,y(n))二0
(1)
2、常微分方程的解析解
有些微分方程可直接通过积分求解.例如,一解常系数常微分方程df=y1可化为
dy=dtt4
y1,两边积分可得通解为y=ce-1.其中c为任意常数.有些常微分方程可用一些技
巧,如分离变量法,积分因子法,常数变异法,降阶法等可化为可积分的方程而求得解析解.
线性常微分方程的解满足叠加原理,从而他们的求解可归结为求一个特解和相应齐次微分方程的通解•一阶变系数线性微分方程总可用这一思路求得显式解。
高阶线性常系数微分方程可用特征根法求得相应齐次微分方程的基本解,再用常数变异法求特解。
一阶常微分方程与高阶微分方程可以互化,已给一个n阶方程
设%=y,y2二y',=y(
y(n)=f(t,y',y"r,y(nJJ)
⑴,可将上式化为一阶方程组
乂'=y2y'n严yn
yn'=f(t,yi,y2,,yn)
反过来,在许多情况下,一阶微分方程组也可化为高阶方程。
所以一阶微分方程组与
高阶常微分方程的理论与方法在许多方面是相通的,一阶常系数线性微分方程组也可用特征根法求解。
3、微分方程的数值解法
除常系数线性微分方程可用特征根法求解,少数特殊方程可用初等积分法求解外,大
部分微分方程无限世界,应用中主要依靠数值解法。
考虑一阶常微分方程初值问题
V(t)=f(t,y(t))t0<^tf
iy(to)=yo(4)
其中y二(yi,y2,…,ym)',f二(fl,f2,…,fm)',y0二(y10,y20,…,ym0)'.所谓数值解法,就是寻求y(t)在一系列离散节点to:
:
:
t^':
:
:
t^-tf上的近似值丫小二0,1,…,n称h/t^-tk为步长,通常取为常量h。
最简单的数值解法是Euler法。
Euler法的思路极其简单:
在节点出用差商近似代替导数
y(tk1)-y(tk)
y'(tk)一-
h
这样导出计算公式(称为Euler格式)
yk1二y「hf(tk,yQ,k二0,1,2,
他能求解各种形式的微分方程。
Euler法也称折线法。
4、解微分方程的MATLAB命令
MATLAB中主要用dsolve求符号解析解,0de45,0de23,ode15s求数值解。
在MATLAB^,由函数dsolve()解决常微分方程(组)的求解问题,其具体格式如下:
r=dsolve('eq1,eq2,...','cond1,cond2,...','v')
'eq1,eq2,...'为微分方程或微分方程组,’cond1,cond2,...',是初始条件或边界条
件,’v'是独立变量,默认的独立变量是't'。
函数dsolve用来解符号常微分方程、方程组,如果没有初始条件,则求出通解,如果有初始条件,则求出特解。
基于龙格-库塔法,MATLAB!
供了求常微分方程的数值解的函数,一般调用格式为:
[t,y]=ode23(filename,tspan,yO)
[t,y]=ode45(filename,tspan,y0)
其中filename是定义f(t,y)的函数文件名,该函数文件必须返回一个列向量。
tspan形
式为[t0,tf],表示求解区间,y0是初始状态列向量。
t和y分别给出时间限量和相应的状
态向量。
ezplot(x,y,[tmin,tmax]):
符号函数的作图命令.x,y为关于参数t的符号函数,[tmin,tmax]为t的取值范围.
inline():
建立一个内联函数.格式:
inline('expr','var1','var2',…),注意括号里的表达式要加引号.
三、实例应用
1、直接用Matlab求微分方程的解析解
22
问题一:
求解微分方程矽/垄的通解.
dx2x2
解:
方程求解的Matlab程序为
y=dsolve('Dy-(xA2+yA2)/xA2/2','x');%求岀的微分方程的解
运行结果为:
y=
x*(-log(x)+2+C1”(-log(x)+C1)
问题二:
求X2鱼2x^ex的通解。
dx
解:
方程求解的Matlab程序为
y=dsolve('Dy*xA2+2*x*y-exp(x)','x')%求岀微分方程的解
运行结果为:
y=
1/xA2*exp(x)+1/xA2*C1
问题三:
求微分方程xy「y-ex=0在初始条件y
(1)=2e下的特解,并画出解函数的图形.
解:
方程求解的Matlab程序为:
symsxy;%定义x,y为符号变量
y=dsolve('x*Dy+y-exp(x)=0','y
(1)=2*exp
(1)','x');%求岀的微分方程在初始条
件y
(1)=2e下的特解
ezplot(y);%作岀解函数的图象
微分方程的特解为:
y=
1/x*exp(x)+1/x*exp
(1)
即yJF,解函解的图形如图1所示。
x
图1解函数的图
2EM
Fiktf匪TgkH■卡
D序口闿A
罟"x-y
问题四:
求q的通解。
先2x—yIdt
解:
方程求解的Matlab程序为:
%求解方程组
[x,y]=dsolve('Dx=4*x-2*y','Dy=2*x-y')
运行结果为:
x=
C1+C2*exp(3*t)y=
1/2*C2*exp(3*t)+2*C1
问题五:
求下列微分方程的解析解
(2)
(3)
方程
(1)求解的MATLAB?
序为:
clear;
s=dsolve('Dy=3*y+4');
运行结果为
清屏
求岀微分方程的解
-4/a+exp(3*t)*C1
方程
(2)求解的MATLAB?
序为:
clear;%
s=dsolve('D2y-sin(2*x)+y','y(0)=0','Dy(0)=1','x')simplify(s)%
运行结果为
s=
清屏
求岀微分方程的解以最简形式显示s
5/3*sin(x)-1/3*sin(2*x)
ans=
-1/3*(-5+2*cos(x))*sin(x)
方程(3)求解的MATLAB?
序为:
clear;%
s=dsolve('Df-f-g','Dg-g+f','f(0)=1','g(0)=1');simplify(s.f);%s
simplify(s.g);
运行结果为:
清屏
求岀微分方程的解
是一个结构
ans=
exp(t)*cos(t)+exp(t)*sin(t)
ans=
-exp(t)*sin(t)+exp(t)*cos(t)
问题六:
利用MATLAB求常微分方程的初值问题(「x2)yJ2xy',y
解。
解:
MATLAB?
序为:
y=dsolve('D2y*(1+xA2)-2*x*Dy=0','y(0)=1,Dy(0)=3','x')
运行结果:
求解微分方程
y=
1+3*x+xA3
⑷
问题七:
利用MATLAB^常微分方程y
-2y'「y''=0的解。
求解微分方程
(组)的初值问题的数值解(近似解)•
解:
MATLAB?
序为:
y=dsolve('D4y-2*D3y+D2y','x')
运行结果为:
y=
C1+C2*x+C3*exp(x)+C4*exp(x)*x
2、用ode23、ode45等求解常微分方程
问题八:
求解微分方程初值问题=_2y2^2x的数值解,求解范围为区间[0,0.5]。
卜(0)=1
解:
MATLAB程序为:
fun=inline('-2*y+2*xA2+2*x','x','y');
%
定义函数
[x,y]=ode23(fun,[0,0.5],1);
%
在区间[0,0.5]
求岀微分方程的数值解
x';
%
输岀y(0)=1
时x的值
y';
%
输岀y(0)=1
时y的值
plot(x,y,'o-');
%
画岀函数的图像
运行结果为:
x'
ans=
0.00000.04000.09000.14000.19000.2400
0.29000.34000.39000.44000.49000.5000y'
ans=
1.00000.92470.84340.77540.71990.6764
0.64400.62220.61050.60840.61540.6179
图形结果如图2所示。
图2函数图象
问题九:
用Euler折线法求解微分方程初值问题
2x
~2,y
的数值解(步长h取0.4),求解范围为区间[0,2].解:
本问题的差分方程为
k=0,1,2,,n-1
(其中:
f(x,y)二y乡).
y
=Xk+h,
y“=yk+hf(Xk,yQI
相应的Matlab程序为:
图3函数图像
问题十:
求解微分方程
y—yt1,y(0)=1,
清屏
求微分方程的解析解以最简形式显示s
先求解析解,再求数值解,并进行比较。
由
clear;%
s=dsolve('Dy=-y+t+1','y(0)=1','t');%
simplify(s)%
运行结果:
ans=
t+exp(-t)
可得解析解为厂t
e
0
下面再求其数值解,先编与M文件fun8.m
M函数fun8.m
functionf=fun8(t,y);
%
编写M文件fun8.m
f=-y+t+1;
%f
的表达式
再用命令
clear;
%
清屏
close;
t=0:
0.1:
1;
%
确定求解区间
y=t+exp(-t);
%
解析解
plot(t,y);
%
画解析解的图形
holdon;
%
保留已经画好的图形,如果下面再画图,两个图形和并在一起
[t,y]=ode45('fun8',[0,1],1);%
调用龙格库塔求解函数求解数值解
plot(t,y);
%
画数值解图形
标岀t,y轴
xlabel('t'),ylabel('y')
结果如图4所示。
图4解析解与数值解
由图4可见,解析解和数值解吻合得很好。
运行结果如图5所示
问题十一:
用欧拉方法和龙格库塔方法求下列微分方程初值问题的数值解,画出解的图形
解:
MATLAB?
序为:
functionf=f(x,y)
%
函数文件
f=y+2*X;
%f
的表达式
clc;clear;
%
清屏
a=0;
%
求解区间
b=1;
[x1,y_r]=ode45('f',[ab],1);
%
调用龙格库塔求解函数求解数值解;
y⑴=1;N=100;h=(b-a)/N;
%
禾U用Euler方法求解
x=a:
h:
b;
%x
的取值
fori=1:
N
%i
的取值
y(i+1)=y(i)+h*f(x(i),y(i));
%y
的表达式
end
%
循环结束
plot(x1,y_r,'r*',x,y,'b+',x,3*exp(x)-2*x-2,'k-');%
画出数值解与真解图
title('数值解与真解图')
;%
标注图形名称
legend('RK4','Euler','
真解');%
用红色*在折线上标出真解
xlabel('x');ylabel('y');
%
标岀x,y轴
0
图5数值解与真解图
%建立函数文件
%y的表达式
%确定t的范围
%t=0时y的值
%求数值解
%求精确解
%绘出图形并比较
问题十二:
设有初值问题:
',y2-1-2
y=4(t1),0汨「
y(0)=2
试求其数值解,并与精确解相比较(精确解为y(t),t11)。
解:
(1)建立函数文件funt.m
functiony=funt(t,y)y=(yA2-t-2)/4/(t+1);
(2)求解微分方程
tO=O;tf=1O;
y0=2;
[t,y]=ode23('funt',[tO,tf],yO);
y1=sqrt(t+1)+1;
plot(t,y,b',t,y1,'r-');
运行结果如图6所示。
图6微分方程精确解与数值解的比较
问题十三:
已知一个二阶线性系统的微分方程为:
x(0)=0,x(0)=1
其中a=2,绘制系统的时间响应曲线和相平面图
解:
令X2=X,Xi二X,则得到系统的状态方程;
X2=X1
I
x1=-ax2
<2(0)=0,X1(0)=1
建立一个函数文件sys.m:
functionxdot=sys(t,x)xdot=[-2*x
(2);x
(1)];
取t0=0,tf=20,求微分方程;
t0=0;tf=20;
[t,x]=ode45('sys',[t0,tf],[1,0]);
[t,x]
subplot(1,2,1);plot(t,x(:
2));subplot(1,2,2);plot(x(:
2),x(:
1));axisequal
%建立函数文件
%xdot的表达式
%确定t的值
%求数值解
%输出结果
%以子图形式绘出解的曲线
%以子图形式绘出相平面曲线
运行结果为:
ans=
0
1.0000
0
0.0001
1.0000
0.0001
0.0001
1.0000
0.0001
0.0002
1.0000
0.0002
0.0002
1.0000
0.0002
0.0005
1.0000
0.0005
0.0007
1.0000
0.0007
0.0010
1.0000
0.0010
0.0012
1.0000
0.0012
19.1332
-0.3498
0.6618
19.2670
-0.5196
0.6036
19.4007
-0.6708
0.5238
19.5344
-0.7980
0.4253
19.6681
-0.8968
0.3116
19.7511
-0.9422
0.2352
19.8340
-0.9747
0.1556
19.9170
-0.9937
0.0738
20.0000
-0.9991
-0.0090
方程的时间相应及相平面曲线如图7所示。
图7时间相应及相平面曲线
四、结论
利用MATLA语言求解微分方程的分析我们不难得出以下结论:
1利用MATLAB勺dsolve函数求微分方程的解析解非常方便快捷。
2、龙格-库塔法求常微分方程的数值解准确方便。
3、在求解常微分方程的过程中,函数文件的使用可是程序简单许多,用起来方便快捷。
4、MATLABt很强大的数学计算和绘图功能,把MATLAB勺的数学计算功能与绘图功能充分的应用,可以方便快捷的解决很多数学方面的问题。
5、在应用MATLAB勺过程中,不但可以普及计算机的利用,充分地利用资源,更可以极大的提高学习者的兴趣,为其他学科的学习提供了较好的计算工具,并且对未来更进一的学习打下了良好的基础。
五、课程体会
1通过运用MATLA可以很轻松的计算出一些方程(组)的解,绘制图形,以及更多的如对高等数学上的应用等。
只要在经过合理的编排后运用MATLAB软件运行即可得出结
果,如果所得结果并非自己想要的,还可以随时修改,直到得出想要结果,如本文的例题,输入程序后,经过检查及修改无误后,再运行程序既可在瞬间的到准确的图形,使用起来也非常方便快捷。
2、通过在本次写论文的过程中,我真的是受益匪浅。
以前也写过不少的论文,但也没有注意论文正规的写作格式,所以在写本篇论文的过程中,我不仅清楚地知道了论文的写作格式,还学会了使用word软件中许多特殊的用法。
在我选择MATLA论文的题目的过程中,我通过查找大量的资料,才发现MATLA涉及的范围非常的大,例如在高等数学中的应用,在大学物理中的应用,在线性代数中的应用等,我才感受到MATLA对我们是如此的
实用。
其次,在刚学习的MATLA的时候,我就觉得它和我们大一第一学期学的C语言有很多相似之处,所以在老师讲的过程中,我就争取好好的听每一节课,每次在实验课上我也认真的做好每一个实验的题目,在程序结果运行出来的时候,我真的特别的兴奋,此时我感觉到了MATLA功能的强大性。
就如本篇论文例举的实例,在平时解决此类问题的时候,非常的麻烦,但通过运用MATLAB军决此类问题,就显得非常的简单了。
但遗憾的就是在老师讲MATLA绘图时,我几乎听的是稀里糊涂,做作业时我也只会编写一些简单的二维和三维的绘图程序,对于图形的精细处理和光照处理等一些复杂一点的操作编程我几乎是一点都不懂。
看到老师绘出的那些经过处理的图形,我十分的羡慕,并且听老师说在我们以后所学到的课程中很多地方都要用到MATLAB所以我一定要好好的学习MATLAB争取让它为我以后的学习中提供更大的帮助。
[参考文献]
[1]刘卫国.MATLA程序设计与应用(第二版)[M].北京:
高等教育出版社,2006.
[2]同济大学数学系.高等数学(第六版)[M].北京:
高等教育出版社,2007.