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 Untitled2 R/ C3 D' p6 b3 H) h: Z
clear all;clc;
) P\" `7 D, G/ N1 H$ B# X$ ^1 a
f=@(t)(2*sin(t)*(t<(4*pi)) + 0);
* G0 P0 f: K# D: @7 x p
g=@(t)(0+cos(t).*(t>=((7*pi)/2)));
3 ^: ~+ V+ H\" l6 h\" {- S- [+ `
function dy = rigid(t,y)7 l; P3 C% \\" \% A6 h) J) Y; \
dy = zeros(2,1);; c* x) ?. \5 @/ C1 A' M- i5 d1 h
dy(1) = y(2)-f(t);
$ g3 e' L( J4 J2 v9 ~; I- u
dy(2) = y(1)*g(t)-y(2);
* D' ^4 E$ [2 { E, G% `
end
* H A3 E, Q+ P) e4 [
options = odeset('RelTol',1e-4,'AbsTol',[1e-5 1e-5]);1 n- {! f2 G7 ?/ y
%[T, y] = ode45(@rigid, [0 20], [1 2],options)9 d7 d\" r\" D. j+ _8 I
sol = ode45(@rigid, [0 20], [1 2],options);
+ k\" m+ ^9 a4 k0 F
x=linspace(0,20,200000);
. D7 X% j; j6 v\" v
y=deval(sol,x);5 h) [7 N/ i* W/ }! n
res=y(1,:)+y(2,:);
; y1 [% t! U+ h$ _\" ?8 q
idx=find(abs(res-0)<1e-4) %相加，当和小于误差运行范围的时候可以认为它就为0
7 G3 Y7 g1 F9 \6 |4 N
xx=x(idx) %算出在x中的下标
. I5 b9 _0 u6 {+ r- l. i1 @
F=@(t)(f(t)+g(t));
2 \5 G2 n0 v- i4 C8 ]
r=[]; %得出的解反代入方程求值，得到的值保存在r中。如果解正确，r中的值应该非常接近0
* E a# ]6 t- P- R& E3 X; @* V. r. v
for i=1:size(xx,2)
; x, x; K8 z9 p5 A2 h, G: s/ C; a
r = [r F(xx(i))]; %将解的值依次带入。当前问题在于r中的值都和0差的很远
9 _) V- j/ F\" i9 f! U( ~. ^( J
end
& o8 x* g; m' B# P\" V
r) e$ m* r( G S* d
end