matlab求解微分方程结果不对 - 数学建模社区-数学中国

function Untitled c& Z% s4 m9 U+ q- `, j9 k! V4 k
clear all;clc;
! y4 {) W/ n$ q" W
f=@(t)(2*sin(t)*(t<(4*pi)) + 0);& C' j5 n+ M4 q1 q
g=@(t)(0+cos(t).*(t>=((7*pi)/2)));8 J9 q; | b) _ J x$ A6 @
function dy = rigid(t,y)
$ x! E$ `& C9 h9 M. L. P5 w& ]) H
dy = zeros(2,1);7 W6 K% k& n8 v; @& }0 L, C
dy(1) = y(2)-f(t);8 m7 [4 U0 u A) P: e3 O! s* B" P$ O
dy(2) = y(1)*g(t)-y(2);1 I" ]0 C# X# F- @* V! k
end
, v% A2 o& g! c" }
options = odeset('RelTol',1e-4,'AbsTol',[1e-5 1e-5]);
: \+ z3 U2 d2 r3 `( Y/ X+ ~% J
%[T, y] = ode45(@rigid, [0 20], [1 2],options). _/ X6 G9 S8 L; _; ]" i
sol = ode45(@rigid, [0 20], [1 2],options);
% |6 j# C; H7 L P& n6 H- O
x=linspace(0,20,200000);
( k P7 l" {+ r( i
y=deval(sol,x);, C; ~, l( ?3 I
res=y(1,:)+y(2,:);6 l2 o' D2 ~* l5 [( B5 I: E
idx=find(abs(res-0)<1e-4) %相加，当和小于误差运行范围的时候可以认为它就为0
+ g# q! ?, O) k+ J+ A, ^. g
xx=x(idx) %算出在x中的下标
2 n v* K, F8 |: j: I7 b6 ?
F=@(t)(f(t)+g(t));
) |+ T% n2 x. a
r=[]; %得出的解反代入方程求值，得到的值保存在r中。如果解正确，r中的值应该非常接近0: S8 Q8 o6 }# h( h
for i=1:size(xx,2)1 g; W. n2 g: {$ T/ A# L2 h" }; Z
r = [r F(xx(i))]; %将解的值依次带入。当前问题在于r中的值都和0差的很远
# M: n; |- g) U {
end
* g( C$ r! O) ?& j2 L4 \: b
r
0 ?- Z) x$ T6 m s% R0 ^9 ?
end

复制代码