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 Untitled5 n' O0 ]; H& i2 e! T
clear all;clc;, R5 Y* H& \- W s8 a, h\" c+ v
f=@(t)(2*sin(t)*(t<(4*pi)) + 0);
: m6 ^# o& Q' A0 u! K9 J
g=@(t)(0+cos(t).*(t>=((7*pi)/2)));
$ x- c7 {; C. E% i }. b
function dy = rigid(t,y)0 i* j3 y& b8 W. ]. m! k
dy = zeros(2,1);4 x K5 b* V+ X+ V0 U- O
dy(1) = y(2)-f(t);5 [: k! e\" G2 n) t9 K) N: S( h
dy(2) = y(1)*g(t)-y(2);
& d& C. M) R; Z- u4 \3 Q! o6 |* p
end
, W% }; k, l, {6 g3 N/ S
options = odeset('RelTol',1e-4,'AbsTol',[1e-5 1e-5]);
0 U! ?& _' G; F& j' A4 g4 K
%[T, y] = ode45(@rigid, [0 20], [1 2],options)\" X& y6 _( R( r |
sol = ode45(@rigid, [0 20], [1 2],options);
/ M* J0 o0 l0 Q/ f/ E
x=linspace(0,20,200000);
0 t1 Z0 \ O) v+ p1 n
y=deval(sol,x);
; X% H5 h9 y5 A
res=y(1,:)+y(2,:);
! \% d+ [& {! i4 @ M
idx=find(abs(res-0)<1e-4) %相加，当和小于误差运行范围的时候可以认为它就为0
$ T0 E' g3 [9 M. j2 |+ S4 K4 k
xx=x(idx) %算出在x中的下标
5 C2 z3 S c& I# |; x
F=@(t)(f(t)+g(t));$ P& i/ `- {6 M, J\" R; i; G6 p
r=[]; %得出的解反代入方程求值，得到的值保存在r中。如果解正确，r中的值应该非常接近0
! m\" R- W. I0 H! r6 ?$ L/ E
for i=1:size(xx,2)) D& n( b; `$ x7 d c9 m3 c
r = [r F(xx(i))]; %将解的值依次带入。当前问题在于r中的值都和0差的很远
2 j3 }9 t5 U* z, C\" L/ }) _
end$ r0 O3 P! L' b# U+ _2 T E* w
r. Q8 D/ e2 w0 P
end