matlab在科学计算中的应用5-3

f (t , x) : 通常微分方程的非齐次项； M(t , x) : 奇异矩阵 为微分代数方程
30
例：
31
• 编写函数
function dx=c5eqdae(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];
& & x1 = x, x2 = x, x3 = y, x4 = y, dx = &&, dy = && x gt; [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-5+x4*x1-x3)/(3*x2+2*x4) dy = (2*x4^2*x1+5-x4*x1+x3)/(3*x2+2*x4)
对于更复杂的问题来说，手工变换的难度将很大，所 以如有可能，可采用计算机去求解有关方程，获得解 析解。如不能得到解析解，也需要在描写一阶常微分 方程组时列写出式子，得出问题的数值解。 17
5.3特殊微分方程的数值解 5.3.1 刚性微分方程的求解
• 刚性微分方程 －刚性过程：微分方程描述的变化过程中，若包含着多 个相互作用但变化速度相差十分悬殊的子过程，这样一 类过程就认为具有 “刚性”。 －刚性微分方程：描述“刚性”过程的微分方程称为 刚性微分方程，相应的初值问题称为“刚性问题” 。 MATLAB采用求解函数ode15s()，该函数的调用格式 和ode45()完全一致。 [t,x]=ode15s(Fun,[t0,tf],x0,options,p1,p2,…)
x1' = x2 ' 2 x2 = − µ ( x1 − 1) x2 − x1
•
>> x0=[2;0]; t_final=3000; >> mu=1000; [t,y]=ode45('vdp_eq',[0,t_final],x0,[],mu); 由于变步长所采用的步长过小，所需时间较长，导致输出的y矩 阵过大，超出计算机存储空间容量。所以不适合采用ode45()来 求解，可用刚性方程求解算法ode15s( )。
18
• 例： ％计算 >> 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 = 1.2970
= y ; x2 = y
'
>> 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),':') 3
34
14
• 若两个高阶微分方程同时含有最高价导数项，需要 先进行相应处理，再应用上述变换方法 • 例：试将二元方程组：
求解： • 先消去其中一个高阶导数
由第一个式子可得：&& = yx + && / 2 y & x
&& = 2 y + 10 − 2 xy − 6 xxy & && x 带入第二个式子得： & & 2 y + 3x
8
• 求解： >> x0=[1.2; 0; 0; -1.04935751];% 输入初值 >> tic, [t,y]=ode45('apolloeq',[0,20],x0); toc elapsed_time = 0.5000 >> length(t), >> plot(y(:,1),y(:,3)) ans = 689
>> [t,x]=ode15s(@c5eqdae,[0,20],x0,options); plot(t,x)
％注意这里使用ode45()函数求解会出错 会
32
•
编写函数： function dx=c5eqdae1(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)];
15
所以：
& x2 = 2 x3 + 10 − 2 x1 x4 − 6 x1 x2 x4
& x = x1 , x = x2 , & y = x3 , y = x4
2 x4 + 3x2
从而得：
&4 = x3 + 5 − x1 x4 + 2 x1 x4 2 x
2 x4 + 3 x2
16
• 用MATLAB符号工具箱求解， ％ 令
23
• 方法二，用ode15s()代替ode45()
>> opt=odeset; opt.RelTol=1e-6; >> tic,[t1,y1]=ode15s('c5exstf2',[0,100],[0;1],opt); toc elapsed_time = 0.2340 >> length(t1), >> plot(t1,y1), legend(‘t1’,’y1’) ans = 169
10
• >> min(diff(t1)) ans = 1.8927e-004 • >> plot(t1(1:end-1),… diff(t1))
11
• 例：
12
• >> x0=[1.2; 0; 0; -1.04935751]; >> tic, [t1,y1]=rk_4('apolloeq',[0,20,0.01],x0); toc elapsed_time = 0.8590 >> plot(y1(:,1),y1(:,3)) % 绘制出轨迹曲线 显而易见，这样求解 是错误的，应该采用 更小的步长。
26
• 应用matlab函数：ode15i
－格式: 格式
%如果不能直接确定初值，可用函数 decic()得出相应初值
27
例：
28
29
5.3.3 微分代数方程求解
• 微分代数方程(differential algebra equation)：指微分 方程中某些变量间满足相应代数方程的约束
• 设
21
• 方法一
>> tic,[t2,y2]=ode45('c5exstf2',[0,100],[0;1]); toc elapsed_time = 33.0630 >> length(t2), plot(t2,y2) ans = 356941
22
• 步长分析：
>> format long, [min(diff(t2)), max(diff(t2))] ans = 0.00022220693884 0.00214971787184 >> plot(t2(1:end-1),diff(t2))
33
>> x0=[0.8; 0.1]; ％也可应用 inline()函数 >> 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')) ％注意这里使 用ode45()函 不会 数求解不会 出错
tocelapsedtime124380计算时间过长绘制出轨迹曲线严格说来某些点仍不满足106的误差限所以求解常微分方程组时建议采用变步长算法而不是定步长算法
5.2.4 微分方程转换
5.2.4.1 单个高阶常微分方程处理方法
1
2
• 例：
• 函数描述为： x1 function y=vdp_eq(t,x,flag,mu) y=[x(2); -mu*(x(1)^2-1)*x(2)-x(1)]; •
24
6.3.2 隐式微分方程求解
• 隐式微分方程为不能转化为显式常微分方程组的方程
例：
25
• 编写函数： function dx=c5ximp(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('c5ximp',[0,10],[0; 0],opt); plot(t,x), 求解过程中Matlab 求解过程中 若提示矩阵奇异 的信息， 的信息，则得出 的数值解无意义. 的数值解无意义
13
• >> tic, [t2,y2]=rk_4('apolloeq',[0,20,0.001],x0); toc elapsed_time = 12.4380 ％计算时间过长 >> plot(y2(:,1),y2(:,3)) % 绘制出轨迹曲线 严格说来某些点仍不 满足10－6的误差限， 所以求解常微分方程 组时建议采用变步长 算法，而不是定步长 算法。
4
6.2.4.2 高阶常微分方程组的变换方法
5
• 例：
6
•
7
• 描述函数： 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];
得出的轨道不正确， 默认精度RelTol设置 得太大，从而导致的 误差传递，可减小该 值。
9
• 改变精度：
>> options=odeset; options.RelTol=1e-6; >> tic, [t1,y1]=ode45('apolloeq',[0,20],x0,options); toc elapsed_time = 0.2970 >> length(t1), >> plot(y1(:,1),y1(:,3)), ans = 1873
19
％作图 >> plot(t,y(:,1)); figure; plot(t,y(:,2)) ％Vandpol方程输出结果y=(y,y’) y(:,1): 解曲线；部分曲线变化较平滑，某些点上变 化在较快, y(:,2) ：解的变化率
20
• 例：
定义函数
function dy=c5exstf2(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];