数学建模社区-数学中国

标题: matlab求解微分方程结果不对 [打印本页]
作者: yzh07137    时间: 2015-7-30 20:35
标题: matlab求解微分方程结果不对
问题是这样的（见附图）9 j) _8 r$ d4 C0 h" _, d

* E+ P8 w* m# b然后我的程序如下
  1. function Untitled
    6 v' E( V5 z2 m' d9 L" ~
  2.         clear all;clc;& ^/ j0 k2 n# C
  3.         f=@(t)(2*sin(t)*(t<(4*pi)) + 0);5 D* U, d# ]6 \) y& Y, Y0 g
  4.         g=@(t)(0+cos(t).*(t>=((7*pi)/2)));8 v: e) U* e; v# r9 {0 a
  5.             function dy = rigid(t,y)$ `/ |) u: j# }% }$ k# Q/ c2 r
  6.                         dy = zeros(2,1);
    & j2 R0 v$ k" h  s
  7.                         dy(1) = y(2)-f(t);
    " i. d4 Q$ K! z/ {! p5 q/ p
  8.                         dy(2) = y(1)*g(t)-y(2);
    : q  g* `: g6 n  R
  9.             end* ]! r% j( o2 D- ?$ Q
  10.             options = odeset('RelTol',1e-4,'AbsTol',[1e-5 1e-5]);2 X" |# R# B, E% l
  11.             %[T, y] = ode45(@rigid, [0 20], [1 2],options)/ z  M, @! W3 A% M/ B
  12.             sol = ode45(@rigid, [0 20], [1 2],options);( }, U3 @( M% x6 g
  13.             x=linspace(0,20,200000);
    , z' g* m! ], u5 x- _# j  S
  14.             y=deval(sol,x);$ A5 ?( H* \2 S; G
  15.             res=y(1,:)+y(2,:);
    0 v8 }! Z, q$ l
  16.             idx=find(abs(res-0)<1e-4) %相加，当和小于误差运行范围的时候可以认为它就为0# @% K/ Y- K/ S! Z
  17.             xx=x(idx) %算出在x中的下标
    5 B0 E( e) ^/ k: e- B% g+ f, q* d
  18.             F=@(t)(f(t)+g(t));) p, g$ v; p9 V1 N) N7 X9 i9 y) @
  19.             r=[]; %得出的解反代入方程求值，得到的值保存在r中。如果解正确，r中的值应该非常接近0
    2 o) \* W, W! N1 j6 @
  20.             for i=1:size(xx,2)
      ?/ g0 w' O* S; `, O  W  E0 [
  21.                     r = [r F(xx(i))]; %将解的值依次带入。当前问题在于r中的值都和0差的很远1 ?$ Z+ v- z1 j5 H+ [2 e( a% i  ]
  22.             end) H& n% g: {' L9 S) W
  23.             r9 Q6 T% X8 i" G9 r8 U  J  K
  24. end
复制代码
问题就是最后算出来的解再带入方程进行检验得出的结果和0差的好远，都到了1.73几。0 ~- P/ c7 p: B" {' m, ]7 v: }
我也不知道哪里用错了，但猜测可能是由于微分方程里出现了f(t)，g(t)的缘故
/ ~1 k* W1 c- c  _  C6 }, F  c$ F* L! @
求教大家我的程序是哪里出现了问题，该如何改呢？8 G8 o9 @" A. U
谢谢- y- n% W8 q# K: i5 O

0 |4 `7 a( t. t- o" `0 [/ M6 \# P, t# d9 S) y! g

matlab.png (12.23 KB, 下载次数: 547)

问题

问题

作者: yzh07137    时间: 2015-8-4 20:45
算了，就知道这种论坛一般也没人
+ G8 J3 s7 p* U* n; m7 p



欢迎光临 数学建模社区-数学中国 (http://www.madio.net/) Powered by Discuz! X2.5