matlab求解微分方程结果不对

字体大小: 正常放大

yzh07137

3 主题	12 听众	150 积分

升级 25%

TA的每日心情

	开心 2015-9-11 00:19

签到天数: 54 天

[LV.5]常住居民I

电梯直达

1^#

发表于 2015-7-30 20:35 |只看该作者 |正序浏览

|招呼Ta 关注Ta

问题是这样的（见附图）

然后我的程序如下

function Untitled
. q\" x\" a1 f2 f$ u# v3 O5 L/ J$ |
clear all;clc;
( z) N1 ]- m, y1 u3 b% D. {& T0 ]
f=@(t)(2*sin(t)*(t<(4*pi)) + 0);4 x+ s: j0 c8 x2 }2 O- v: p9 W
g=@(t)(0+cos(t).*(t>=((7*pi)/2)));
% R. o! N) \( ^3 S& w\" I
function dy = rigid(t,y)
: o* D7 A; V\" [+ \) e
dy = zeros(2,1);
- H3 e. M* Q9 Z# y) n) z8 q
dy(1) = y(2)-f(t);
( l) y) Y. O8 Q\" C! M! i: e; {. w8 K
dy(2) = y(1)*g(t)-y(2);0 i; V* Y8 v7 {
end: r: W$ h. u% @1 \4 R
options = odeset('RelTol',1e-4,'AbsTol',[1e-5 1e-5]);/ y7 M& s+ ~3 U8 i8 ?* o
%[T, y] = ode45(@rigid, [0 20], [1 2],options)
9 t: D: F5 ^/ C0 y5 x( x* P
sol = ode45(@rigid, [0 20], [1 2],options);
5 r% j( a6 n) U\" G$ F
x=linspace(0,20,200000);; r8 S% D2 _ u; W+ b. ^
y=deval(sol,x);) z. @: }4 V, G
res=y(1,:)+y(2,:);
8 A {2 z! [9 E% G
idx=find(abs(res-0)<1e-4) %相加，当和小于误差运行范围的时候可以认为它就为0; }, M( R: s- I5 _+ \- J1 y3 j' e
xx=x(idx) %算出在x中的下标
+ \/ B7 R6 T/ r
F=@(t)(f(t)+g(t));
1 o! h' ^) l; h. E. n2 H' K4 I- X
r=[]; %得出的解反代入方程求值，得到的值保存在r中。如果解正确，r中的值应该非常接近0' K* L2 \- \' O: g' R0 s
for i=1:size(xx,2)
) s2 \& |0 l4 [( r( c
r = [r F(xx(i))]; %将解的值依次带入。当前问题在于r中的值都和0差的很远. k4 ?% \4 G2 @/ l$ @\" E
end, S/ d6 v% {9 n' l& ~1 y g: A
r
6 H# l\" a9 ?4 \; ]
end