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
: s; W/ ?7 g/ h6 f) E
clear all;clc;8 D) B. N4 g. K x
f=@(t)(2*sin(t)*(t<(4*pi)) + 0);- T) [5 C) H5 [' N. J+ Q! G( G+ V
g=@(t)(0+cos(t).*(t>=((7*pi)/2)));
4 s4 s4 E: }, ^
function dy = rigid(t,y) g4 Q& U1 B$ o. z/ p: S* O# b
dy = zeros(2,1);$ w! r; T) O' e8 a5 T3 o3 g9 Q
dy(1) = y(2)-f(t);
! ~! G) _1 Q+ P& _
dy(2) = y(1)*g(t)-y(2);
9 T2 z; A- {9 S) W& a7 D& y2 K
end
; g' x- A: Q f- b c7 f! B8 J' y
options = odeset('RelTol',1e-4,'AbsTol',[1e-5 1e-5]);) ^7 w* t2 Q6 I/ d* b8 Q! v
%[T, y] = ode45(@rigid, [0 20], [1 2],options)
# ~* o0 X+ ?$ |
sol = ode45(@rigid, [0 20], [1 2],options);, m: U, {+ S6 h3 ]* s
x=linspace(0,20,200000);
3 w% v9 u, G u$ @; b* c
y=deval(sol,x);
# k! s6 u2 ^+ Q) ^( n+ w; J
res=y(1,:)+y(2,:);6 f/ D\" K% m0 F+ i3 ?1 ?5 _
idx=find(abs(res-0)<1e-4) %相加，当和小于误差运行范围的时候可以认为它就为0
1 P P$ t: x, ]- q
xx=x(idx) %算出在x中的下标/ ~7 e3 y. ?( ?: c8 O# [
F=@(t)(f(t)+g(t));
8 R0 l- s9 w/ S, a( y4 C, {, I' V
r=[]; %得出的解反代入方程求值，得到的值保存在r中。如果解正确，r中的值应该非常接近0
! P6 u9 Y/ e/ h1 W) B7 V
for i=1:size(xx,2)3 a3 F! E( x' |4 Q0 d- ?2 R5 [
r = [r F(xx(i))]; %将解的值依次带入。当前问题在于r中的值都和0差的很远7 h8 o6 D8 b$ T6 i0 Z
end# \8 `9 |+ q' l2 }; L& ~- p
r\" A; |* K% \6 p; ~
end