koorio.com
海量文库 文档专家
当前位置:首页 >> 数学 >>

MALAB微分方程问题的解法_图文

第六章 微分方程问题的解法
? 微分方程的解析解方法 ? 常微分方程问题的数值解法
–微分方程问题算法概述 –四阶定步长 Runge-Kutta算法及 MATLAB 实现 –一阶微分方程组的数值解 –微分方程转换

? 特殊微分方程的数值解 ? 边值问题的计算机求解 ? 偏微分方程的解

6.1 微分方程的解析解方法
? 格式: y=dsolve(f1, f2, …, fm) ? 格式:指明自变量 y=dsolve(f1, f2, …, fm ,’x’) fi即可以描述微分方程,又可描述初始条件 或边界条件。如: 描述微分方程时 y
(4)

(t ) ? 7 ? D4 y ? 7

描述条件时

?? y(2) ? 3 ? D2 y(2) ? 3

例:

>> syms t; u=exp(-5*t)*cos(2*t+1)+5; >> uu=5*diff(u,t,2)+4*diff(u,t)+2*u uu = 87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+10 >> syms t y; >> y=dsolve(['D4y+10*D3y+35*D2y+50*Dy+24*y=',... '87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+10'])

>> y=dsolve(['D4y+10*D3y+35*D2y+50*Dy+24*y=',... '87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1) ... +10'], 'y(0)=3', 'Dy(0)=2', 'D2y(0)=0', 'D3y(0)=0')

分别处理系数,如: >> [n,d]=rat(double(vpa(-445/26*cos(1)-51/13*sin(1)-69/2)))] ans = -8704 185 % rat()最接近有理数的分数

判断误差: >> vpa(-445/26*cos(sym(1))-51/13*sin(1)-69/2+8704/185) ans = .114731975864790922564144636e-4

>> y=dsolve(['D4y+10*D3y+35*D2y+50*Dy+24*y=',...

'87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1) + ... 10'],'y(0)=1/2','Dy(pi)=1','D2y(2*pi)=0','Dy(2*pi)=1/5');
如果用推导的方法求Ci的值,每个系数的解析解至少要写出10 数行,故可采用有理式近似 的方式表示. >> vpa(y,10) %有理近似值 ans = 1.196361839*exp(-5.*t)+.4166666667.4785447354*sin(t)*cos(t)*exp(-5.*t)-.4519262218e1*cos(2.*t)*exp(-5.*t)-2.392723677*cos(t)^2*exp(5.*t)+.2259631109*sin(2.*t)*exp(-5.*t)-473690.0893*exp(3.*t)+31319.63786*exp(-2.*t)-219.1293619*exp(1.*t)+442590.9059*exp(-4.*t)

? 例:求解

>> [x,y]=dsolve('D2x+2*Dx=x+2*y-exp(-t)', … 'Dy=4*x+3*y+4*exp(-t)')

? 例:
>> syms t x >> x=dsolve('Dx=x*(1-x^2)') x= [ 1/(1+exp(-2*t)*C1)^(1/2)] [ -1/(1+exp(-2*t)*C1)^(1/2)]

>> syms t x; x=dsolve('Dx=x*(1-x^2)+1') Warning: Explicit solution could not be found; implicit solution returned. > In D:\MATLAB6p5\toolbox\symbolic\dsolve.m at line 292 x= t-Int(1/(a-a^3+1),a=``..x)+C1=0 故只有部分非线性微分方程有解析解。

6.2 微分方程问题的数值解法
6.2.1 微分方程问题算法概述

微分方程求解的误差与步长问题:

6.2.2 四阶定步长Runge-Kutta算法 及 MATLAB 实现

function [tout,yout]=rk_4(odefile,tspan,y0) %y0初值列向量 t0=tspan(1); th=tspan(2); if length(tspan)<=3, h=tspan(3); % tspan=[t0,th,h] else, h=tspan(2)-tspan(1); th=tspan(end); end %等间距数组 tout=[t0:h:th]'; yout=[]; for t=tout' k1=h*eval([odefile ‘(t,y0)’]); % odefile是一个字符串变 量,为表示微分方程f( )的文件名。 k2=h*eval([odefile '(t+h/2,y0+0.5*k1)']); k3=h*eval([odefile '(t+h/2,y0+0.5*k2)']); k4=h*eval([odefile '(t+h,y0+k3)']); y0=y0+(k1+2*k2+2*k3+k4)/6; yout=[yout; y0']; end %由效果看,该算法不是一个较好的方法。

6.2.3 一阶微分方程组的数值解
6.2.3.1 四阶五级Runge-Kutta-Felhberg算法

通过误差向量调节步长,此为自动变步长方法。 四阶五级RKF算法有参量系数表。

6.2.3.2 基于 MATLAB 的微分方程
求解函数 格式1: 直接求解 [t,x]=ode45(Fun,[t0,tf],x0)

格式2: 带有控制参数 [t,x]=ode45(Fun,[t0,tf],x0,options) 格式3: 带有附加参数 [t,x]=ode45(Fun,[t0,tf],x0,options,p1,p2,…)

[t0,tf]求解区间, x0初值问题的初始状态变量。

描述需要求解的微分方程组: 不需附加变量的格式 function xd=funname(t,x) 可以使用附加变量 function xd=funname(t,x,flag,p1,p2,…) % t是时间变量或自变量(必须给),x为状态向量, xd为返回状态向量的导数。flag用来控制求解过程,指
定初值,即使初值不用指定,也必须有该变量占位。

修改变量:options唯一结构体变量,用odeset( )修改。 options=odeset(‘RelTol’,1e-7); options= odeset; options. RelTol= 1e-7;

? 例:

自变函数
function xdot = lorenzeq(t,x) xdot=[-8/3*x(1)+x(2)*x(3); -10*x(2)+10*x(3);… -x(1)*x(2)+28*x(2)-x(3)];

>> t_final=100; x0=[0;0;1e-10]; % t_final为设定的仿真终止时间 >> [t,x]=ode45('lorenzeq',[0,t_final],x0); plot(t,x), >> figure; % 打开新图形窗口 >> plot3(x(:,1),x(:,2),x(:,3)); >> axis([10 42 -20 20 -20 25]); % 根据实际数值手动设 置坐标系

? 可采用comet3( )函数绘制动画式的轨迹。 >> comet3(x(:,1),x(:,2),x(:,3))

? 描述微分方程是常微分方程初值问题数值求解 的关键。
>> f1=inline(['[-8/3*x(1)+x(2)*x(3); -10*x(2)+10*x(3);',... '-x(1)*x(2)+28*x(2)-x(3)]'],'t','x'); >> t_final=100; x0=[0;0;1e-10]; >> [t,x]=ode45(f1,[0,t_final],x0); >> plot(t,x), figure; >> plot3(x(:,1),x(:,2),x(:,3)); axis([10 42 -20 20 -20 25]); 得出完全一致的结果。

6.2.3.3 MATLAB 下带有附加参数的微 分方程求解
? 例:

? 编写函数 function xdot=lorenz1(t,x,flag,beta,rho,sigma)
% flag变量是不能省略的

xdot=[-beta*x(1)+x(2)*x(3); -rho*x(2)+rho*x(3); -x(1)*x(2)+sigma*x(2)-x(3)]; 求微分方程: >> t_final=100; x0=[0;0;1e-10]; >> b2=2; r2=5; s2=20; >> [t2,x2]=ode45('lorenz1',[0,t_final],x0,[],b2,r2,s2); >> plot(t2,x2), %options位置为[],表示不需修改控制选项 >> figure; plot3(x2(:,1),x2(:,2),x2(:,3)); axis([0 72 -20 22 -35 40]);

f2=inline(['[-beta*x(1)+x(2)*x(3); -rho*x(2)+rho*x(3);',... '-x(1)*x(2)+sigma*x(2)-x(3)]'], … 't','x','flag','beta','rho','sigma'); % flag变量是不能省略的

6.2.4 微分方程转换
6.2.4.1 单个高阶常微分方程处理方法

? 例:
? 函数描述为: function y=vdp_eq(t,x,flag,mu) y=[x(2); -mu*(x(1)^2-1)*x(2)-x(1)];

>> x0=[-0.2; -0.7]; t_final=20; >> mu=1; [t1,y1]=ode45('vdp_eq',[0,t_final],x0,[],mu); >> mu=2; [t2,y2]=ode45('vdp_eq',[0,t_final],x0,[],mu); >> plot(t1,y1,t2,y2,':') >> figure; plot(y1(:,1),y1(:,2),y2(:,1),y2(:,2),':')

>> x0=[2;0]; t_final=3000; >> mu=1000; [t,y]=ode45('vdp_eq',[0,t_final],x0,[],mu);

由于变步长所采用的步长过小,所需时间较长,导致输出的y 矩阵过大,超出计算机存储空间容量。所以不适合采用 ode45()来求解,可用刚性方程求解算法ode15s( )。

6.2.4.2 高阶常微分方程组的变换方法

? 例:

? 描述函数:
function dx=apolloeq(t,x) mu=1/82.45; mu1=1-mu; r1=sqrt((x(1)+mu)^2+x(3)^2); r2=sqrt((x(1)-mu1)^2+x(3)^2); dx=[x(2); 2*x(4)+x(1)-mu1*(x(1)+mu)/r1^3-mu*(x(1)-mu1)/r2^3; x(4); -2*x(2)+x(3)-mu1*x(3)/r1^3-mu*x(3)/r2^3];

? 求解: >> x0=[1.2; 0; 0; -1.04935751]; >> tic, [t,y]=ode45('apolloeq',[0,20],x0); toc elapsed_time = 0.8310 >> length(t), >> plot(y(:,1),y(:,3)) ans = 689
得出的轨道不正确, 默认精度RelTol设置 得太大,从而导致的 误差传递,可减小该 值。

? 改变精度:
>> options=odeset; options.RelTol=1e-6; >> tic, [t1,y1]=ode45('apolloeq',[0,20],x0,options); toc elapsed_time = 0.8110 >> length(t1), >> plot(y1(:,1),y1(:,3)), ans = 1873

>> min(diff(t1)) ans = 1.8927e-004

>> plot(t1(1:end-1),… diff(t1))

? 例:

>> x0=[1.2; 0; 0; -1.04935751]; >> tic, [t1,y1]=rk_4('apolloeq',[0,20,0.01],x0); toc elapsed_time = 4.2570 >> plot(y1(:,1),y1(:,3)) % 绘制出轨迹曲线

显而易见,这样求解 是错误的,应该采用 更小的步长。

>> tic, [t2,y2]=rk_4('apolloeq',[0,20,0.001],x0); toc elapsed_time = 124.4990 %计算时间过长 >> plot(y2(:,1),y2(:,3)) % 绘制出轨迹曲线
严格说来某些点仍不 满足10-6的误差限, 所以求解常微分方程 组时建议采用变步长 算法,而不是定步长 算法。

? 例:

?, x3 ? y, x4 ? y ? , dx ? ?? x1 ? x, x2 ? x x, dy ? ?? y

用MATLAB符号工具箱求解, 令



>> syms x1 x2 x3 x4 >> [dx,dy]=solve(‘dx+2*x4*x1=2*dy’, ‘dx*x4+ … 3*x2*dy+x1*x4-x3=5’,‘dx,dy’) % dx,dy为指定变量 dx = -2*(3*x4*x1*x2+x4*x1-x3-5)/(2*x4+3*x2) dy = (2*x4^2*x1-x4*x1+x3+5)/(2*x4+3*x2) 对于更复杂的问题来说,手工变换的难度将很 大,所以如有可能,可采用计算机去求解有关方程, 获得解析解。如不能得到解析解,也需要在描写一 阶常微分方程组时列写出式子,得出问题的数值解。

6.3特殊微分方程的数值解 6.3.1 刚性微分方程的求解
? 刚性微分方程 一类特殊的常微分方程,其中一些解变 化缓慢,另一些变化快,且相差悬殊,这类 方程常常称为刚性方程。

MATLAB采用求解函数ode15s(),该函数的调用 格式和ode45()完全一致。
[t,x]=ode15s(Fun,[t0,tf],x0,options,p1,p2,…)

? 例:
%计算 >> h_opt=odeset; h_opt.RelTol=1e-6; >> x0=[2;0]; t_final=3000; >> tic, mu=1000; [t,y]=ode15s('vdp_eq',[0,t_final],x0,h_opt,mu); toc elapsed_time = 2.5240

%作图 >> plot(t,y(:,1)); figure; plot(t,y(:,2))
y(:,1)曲线变化较平滑, y(:,2)变化在某些点上 较快。

? 例:

定义函数
function dy=c7exstf2(t,y) dy=[0.04*(1-y(1))-(1-y(2))*y(1)+0.0001*(1-y(2))^2; -10^4*y(1)+3000*(1-y(2))^2];

? 方法一
>> tic,[t2,y2]=ode45('c7exstf2',[0,100],[0;1]); toc elapsed_time = 229.4700 >> length(t2), plot(t2,y2) ans = 356941

? 步长分析:
>> format long, [min(diff(t2)), max(diff(t2))] ans = 0.00022220693884 0.00214971787184 >> plot(t2(1:end-1),diff(t2))

? 方法二,用ode15s()代替ode45()
>> opt=odeset; opt.RelTol=1e-6; >> tic,[t1,y1]=ode15s('c7exstf2',[0,100],[0;1],opt); toc elapsed_time = 0.49100000000000 >> length(t1), >> plot(t1,y1) ans = 169

6.3.2 隐式微分方程求解
? 隐式微分方程为不能转化为显式常微分方程组的方程

例:

? 编写函数: function dx=c7ximp(t,x) A=[sin(x(1)) cos(x(2)); -cos(x(2)) sin(x(1))]; B=[1-x(1); -x(2)]; dx=inv(A)*B; 求解: >> opt=odeset; opt.RelTol=1e-6; >> [t,x]=ode45('c7ximp',[0,10],[0; 0],opt); plot(t,x)

6.3.3 微分代数方程求解
例:

? 编写函数
function dx=c7eqdae(t,x) dx=[-0.2*x(1)+x(2)*x(3)+0.3*x(1)*x(2); 2*x(1)*x(2)-5*x(2)*x(3)-2*x(2)*x(2); x(1)+x(2)+x(3)-1]; >> M=[1,0,0; 0,1,0; 0,0,0]; >> options=odeset; >> options.Mass=M; % Mass微分代数方程中 的质量矩阵(控制参数) >> x0=[0.8; 0.1; 0.1];

>> [t,x]=ode15s(@c7eqdae,[0,20],x0,options); plot(t,x)

编写函数: function dx=c7eqdae1(t,x) dx=[-0.2*x(1)+x(2)*(1-x(1)-x(2))+0.3*x(1)*x(2); 2*x(1)*x(2)-5*x(2)*(1-x(1)-x(2))-2*x(2)*x(2)];

>> x0=[0.8; 0.1]; >> fDae=inline(['[-0.2*x(1)+x(2)*(1-x(1)x(2))+0.3*x(1)*x(2);',... '2*x(1)*x(2)-5*x(2)*(1-x(1)-x(2))-2*x(2)*x(2)]'],'t','x'); >> [t1,x1]=ode45(fDae,[0,20],x0); plot(t1,x1,t1,1-sum(x1'))

6.3.3延迟微分方程求解

sol:结构体数据,sol.x:时间向量t, sol.y:状态向量。
f1 : 延迟微分方程,? ? [?1 , ?? n ], f 2 : t ? t0时的状态变量值函数。

? 例:

编写函数: function dx=c7exdde(t,x,z) xlag1=z(:,1); %第一列表示提取 x(? 1 ) xlag2=z(:,2); dx=[1-3*x(1)-xlag1(2)-0.2*xlag2(1)^3-xlag2(1); x(3); 4*x(1)-2*x(2)-3*x(3)]; 历史数据函数: function S=c7exhist(t) S=zeros(3,1);

? 求解: >> lags=[1 0.5]; tx=dde23('c7exdde',lags,zeros(3,1),[0,10]); >> plot(tx.x,tx.y(2,:)) %与ode45()等返回的x矩阵不一样,它是按行排列的。

6.4边值问题的计算机求解

6.4.1 边值问题的打靶算法
数学方法描述: 以二阶方程为例

? y '' ? F ( x, y, y ' ) ? ? y(a) ? ? , y(b) ? ?
m1 m2

?1
?

? y '' ? F ( x, y, y ' ) ? ? ' ? y (a) ? ? , y (a) ? m
m2 ? m1 m3 ? m1 ? ( ? ? ?1 ) ? 2 ? ?1

m

?2

F ( x, y, y ) 线性的 编写函数: function [t,y]=shooting(f1,f2,tspan,x0f,varargin) t0=tspan(1); tfinal=tspan(2); ga=x0f(1); gb=x0f(2); [t,y1]=ode45(f1,tspan,[1;0],varargin); [t,y2]=ode45(f1,tspan,[0;1],varargin); [t,yp]=ode45(f2,tspan,[0;0],varargin); m=(gb-ga*y1(end,1)-yp(end,1))/y2(end,1); [t,y]=ode45(f2,tspan,[ga;m],varargin);

'

例:

编写函数: function xdot=c7fun1(t,x) xdot=[x(2); -2*x(1)+3*x(2)];
function xdot=c7fun2(t,x) xdot=[x(2); t-2*x(1)+3*x(2)]; >> [t,y]=shooting('c7fun1', … 'c7fun2',[0,1],[1;2]); plot(t,y)

原方程的解析解为

解的检验 >> y0=((exp(2)-3)*exp(t)+(3-exp(1))*exp(2*t))/(4*exp(1)*(exp(1)1))+3/4+t/2; >> norm(y(:,1)-y0) % 整个解函数检验 ans = 4.4790e-008 >> norm(y(end,1)-2) % 终点条件检验 ans = 2.2620e-008

非线性方程边值问题的打靶算法:

? y ? F ( x, y , y ) ? ? y(a) ? ? , y(b) ? ?
'' '

? y ? F ( x, y , y ) ? ? ' ? y (a) ? ? , y (a) ? m
'' '

用Newton迭代法处理

y '' ? F ( x, y, y ' )

编写函数: function [t,y]=nlbound(funcn,funcv,tspan,x0f,tol,varargin) t0=tspan(1);tfinal=tspan(2); ga=x0f(1); gb=x0f(2); m=1; m0=0; while (norm(m-m0)>tol), m0=m; [t,v]=ode45(funcv,tspan,[ga;m;0;1],varargin); m=m0-(v(end,1)-gb)/(v(end,3)); end [t,y]=ode45(funcn,tspan,[ga;m],varargin);

例:

编写两个函数: function xdot=c7fun3(t,x) xdot=[x(2); 2*x(1)*x(2); x(4); 2*x(2)*x(3)+2*x(1)*x(4)]; function xdot=c7fun4(t,x) xdot=[x(2); 2*x(1)*x(2)];

>> [t,y]=nlbound('c7fun4','c7fun3',[0,pi/2],[-1,1],1e-8); >> plot(t,y); set(gca,'xlim',[0,pi/2]); 精确解:
检验: >> y0=tan(t-pi/4); >> norm(y(:,1)-y0) ans = 1.6629e-005 >> norm(y(end,1)-1) ans = 5.2815e-006

6.4.2 线性微分方程的有限差分算法
y(a) ? ? a , y(b) ? ? b
把等式左边用差商表示。

编写函数: function [x,y]=fdiff(funcs,tspan,x0f,n) t0=tspan(1);tfinal=tspan(2); ga=x0f(1); gb=x0f(2); h=(tfinal-t0)/n; for i=1:n, x(i)=t0+h*(i-1); end, x0=x(1:n-1); t=-2+h^2*feval(funcs,x0,2); tmp=feval(funcs,x0,1); v=1+h*tmp/2; w=1-h*tmp/2; b=h^2*feval(funcs,x0,3); b(1)=b(1)-w(1)*ga; b(n-1)=b(n-1)-v(n-1)*gb; b=b'; A=diag(t); for i=1:n-2, A(i,i+1)=v(i); A(i+1,i)=w(i+1); end y=inv(A)*b; x=[x tfinal]; y=[ga; y; gb]';

例:
编写函数: function y=c7fun5(x,key) switch key case 1, y=1+x; case 2, y=1-x; otherwise, y=1+x.^2; end >> [t,y]=fdiff('c7fun5',[0,1],[1,4],50); plot(t,y)

6.5 偏微分方程求解入门
6.5.1 偏微分方程组求解

函数描述:

? 边界条件的函数描述:

? 初值条件的函数描述:
u0=pdeic(x)

? 例:

? 函数描述:

function [c,f,s]=c7mpde(x,t,u,du) c=[1;1]; y=u(1)-u(2); F=exp(5.73*y)-exp(-11.46*y); s=F*[-1; 1]; f=[0.024*du(1); 0.17*du(2)];

描述边界条件的函数 function [pa,qa,pb,qb]=c7mpbc(xa,ua,xb,ub,t) pa=[0; ua(2)]; qa=[1;0]; pb=[ub(1)-1; 0]; qb=[0;1];

? 描述初值:
function u0=c7mpic(x) u0=[1; 0];

求解: >> x=0:.05:1; t=0:0.05:2; m=0; >> sol=pdepe(m,@c7mpde,@c7mpic,@c7mpbc,x,t); >> surf(x,t,sol(:,:,1)), figure; surf(x,t,sol(:,:,2))

6.5.2 二阶偏微分方程的数学描述
? 椭圆型偏微分方程:

? 抛物线型偏微分方程:

? 双曲型偏微分方程:

? 特征值型偏微分方程:

6.5.3 偏微分方程的求解界面应用简介 6.5.3.1 偏微分方程求解程序概述
? 启动偏微分方程求解界面
– 在 MATLAB 下键入 pdetool

? 该界面分为四个部分
– – – – 菜单系统 工具栏 集合编辑 求解区域

?

6.5.3.2 偏微分方程求解区域绘制
1)用工具栏中的椭圆、矩形等绘制一些区域。 2)在集合编辑栏中修改其内容。 如(R1+E1+E2)-E3 3)单击工具栏中 按纽可得求解边界。 4)选择Boundary-Remove All Subdomain Borders菜单项,消除相邻区域中间的分隔线。 5)单击 ? 按纽可将求解区域用三角形划分成网 格。可用 ? 按纽加密。

??

6.5.3.3 偏微分方程边界条件描述

选择Boundary-Specify Boundary Conditions菜单

6.5.3.4 偏微分方程求解举例
? 例:

求解: 1)绘制求解区域。 2)描述边界条件(Boundary-Specify Boundary Conditions)。 3)选择偏微分方程的类型。单击工具栏中的PDE图标, 在打开的新窗口选择Hyperbolic选项,输入参数c,a,f,d. 4)求解。单击工具栏中的等号按钮。

显示: 1)图形颜色表示t=0时u(x,y)的函数值。 2)单击工具栏中的三维图标将打开一新的对话框,若 再选择Contour可绘制等值线,若选择Arrows选项可 绘制引力线。若单独选择Height(3d-plot),则在另 一窗口绘制出三维图形。 3)可在单击三维图标打开的新对话框中,对Property 栏目的各个项目重新选择。 4)可修改微分方程的边界条件,重新求解。 动画: 1)Solve-Parameters对话框时间向量改为0:0.1:2。 2)三维图标打开的对话框中选择Animation选项,单 击Options按纽设置播放速度。Plot-Export Movie 菜 单。

6.5.3.5 函数参数的偏微分方程求解
? 例:(椭圆型)

? 求解: 1)求解区域不变。 2)描述边界条件,u=0。 3)选择偏微分方程的类型。单击工具栏中的 PDE图标,在打开的新窗口选择Elliptic选项, 输入参数c=1./sqrt(1+ux.^2+uy.^2), a=x.^2+y.^2 , f=exp(-x.^2-y.^2). 4)再打开Solve-Parameters对话框,选定Use nonlinear solve属性(该属性只适于椭圆性偏微 分方程) 5)求解。单击工具栏中的等号按钮。


网站首页 | 网站地图
All rights reserved Powered by 酷我资料网 koorio.com
copyright ©right 2014-2019。
文档资料库内容来自网络,如有侵犯请联系客服。3088529994@qq.com