- 在线时间
- 0 小时
- 最后登录
- 2009-8-17
- 注册时间
- 2009-8-5
- 听众数
- 4
- 收听数
- 0
- 能力
- 0 分
- 体力
- 74 点
- 威望
- 0 点
- 阅读权限
- 20
- 积分
- 43
- 相册
- 1
- 日志
- 2
- 记录
- 2
- 帖子
- 17
- 主题
- 4
- 精华
- 1
- 分享
- 1
- 好友
- 6
群组: 数学建模 |
Matlab中微分方程求解
前面已经介绍了微分方程(组)的求解方法:解析解法、数值解法和图解法等等。这些方法计算工作量大,需要借助于计算机实现。下面简单介绍Matlab在此领域的应用。
3.1 解析解
desolve(‘eqn1’,’eqn2’, ...) 求常微分方程(组)的解析解。
1. 微分方程
例1:求解二阶微分方程x2y’’ + xy’ + (x2-n2)y = 0, y (p/2) = 2, y’(p /2) = - 2/p, n = 1/2.
解析解:
dsolve('D2y+(1/x)*Dy+(1-(1/2)^2/x^2)*y=0','y(pi/2)=2,Dy(pi/2)=-2/pi','x')
ans = 2^(1/2)*pi^(1/2)*sin(x)/x^(1/2)
pretty(ans)
2.微分方程组
例2:求解 df/dx=3f+4g; dg/dx=-4f+3g。
通解: [f,g]=dsolve('Df=3*f+4*g','Dg=-4*f+3*g')
f= exp(3*t)*cos(4*t)*C1+exp(3*t)*sin(4*t)*C2
g=-exp(3*t)*sin(4*t)*C1+exp(3*t)*cos(4*t)*C2
特解: [f,g]=dsolve('Df=3*f+4*g','Dg=4*f+3*g','f(0)=0,g(0)=1')
f=exp(3*t)*sin(4*t)
g=exp(3*t)*cos(4*t)
3.2 数值解
在微分方程(组)难以获得解析解的情况下,可以用Matlab方便地求出数值解。格式为:
[t,y] = ode23('F',TSPAN,Y0,OPTIONS,P1,P2,...)
注意:微分方程组形式:y' = F(t, y),t为自变量,y为因变量(可以是多个,如微分方程组)。[t, y]为输出矩阵,分别表示自变量和因变量的取值;F代表微分方程组的函数名(m文件,必须返回一个列向量);TSPAN=[t0, tf] 表示自变量的取值范围;Y0为初值条件;OPTIONS为选项;P1,P2,...用于传递附加的参数值。
ode23是微分方程组数值解的低阶方法,ode45为中阶方法,与ode23类似。另外还有系列其它方法,如求解非刚性微分方程组的可变阶方法 ode123。
例 3 求解一个经典的(Van Der pol)微分方程:
解 形式转化:令y1=u(t); y2= (t)。则以上方程转化为: 。
编写M文件如下(必须是M文件表示微分方程组,一般地,M文件的名字与函数名相同):
function dot1=vdpol(t,y);
dot1=[y(2); (1-y(1)^2)*y(2)-y(1)];
在工作空间写如下命令:
[t,y]=ode23('vdpol',0,20,[1,0]);y1=y(:,1);y2=y(:,2);
plot(t,y1,t,y2,'--');title('Van Der Pol Solution ');
xlabel('Time,Second');ylabel('y(1)andy(2)')
注:Van der Pol方程描述具有一个非线性振动项的振动子的运动过程。最初,由于它在非线性电路上的应用而引起广泛兴趣。一般形式为 。
又如例2的数值解:
首先建立M-文件函数:
function f=jie3(x,y)
f=[y(1) -y(2)/x+((1/2)^2/x^2-1)*y(1)];
计算:
[x,y]=ode23('jie3',[pi/2,pi],[2,-2/pi])
x = 1.5708 1.6074 1.7645 1.9215 2.0786 2.2357 2.3928 2.5499 2.7069
2.8640 3.0211 3.1416;
y =
2.0000 -0.6366
2.0745 -0.6885
2.4273 -0.9351
2.8401 -1.2259
3.3231 -1.5695
3.8883 -1.9758
4.5495 -2.4564
5.3232 -3.0249
6.2285 -3.6973
7.2877 -4.4923
8.5271 -5.4319
9.6189 -6.2668
作图如下:
y1=y(:,1);
y2=y(:,2);
plot(x,y1,'--',x,y2,'r')
plot(x,y1,x,y2,'r') 图3-8 例2的数值解图形
3.3 图形解
无论是解析解还是数值解,都不如图形解直观明了。即使是在得到了解析解或数值解的情况下,作出解的图形,仍然是一件深受欢迎的事。这些都可以用Matlab方便地进行。
1.图示解析解
如果微分方程(组)的解析解为:y=f (x),则可以用Matlab函数fplot作出其图形:
fplot(fun,lims)
其中:fun给出函数表达式;lims=[xmin xmax ymin ymax]限定坐标轴的大小。例如
fplot('sin(1/x)', [0.01 0.1 -1 1])
2.图示数值解
设想已经得到微分方程(组)的数值解(x,y)。可以用Matlab函数plot(x,y)直接作出图形。其中x和y为向量(或矩阵)。 |
|