数学建模社区-数学中国

标题: matlab求解微分方程结果不对 [打印本页]

作者: yzh07137    时间: 2015-7-30 20:35
标题: matlab求解微分方程结果不对
问题是这样的(见附图)
* B, ?7 s  V) |! J* a
3 H% t* w6 f& h! A; P然后我的程序如下
  1. function Untitled  c& Z% s4 m9 U+ q- `, j9 k! V4 k
  2.         clear all;clc;
    ! y4 {) W/ n$ q" W
  3.         f=@(t)(2*sin(t)*(t<(4*pi)) + 0);& C' j5 n+ M4 q1 q
  4.         g=@(t)(0+cos(t).*(t>=((7*pi)/2)));8 J9 q; |  b) _  J  x$ A6 @
  5.             function dy = rigid(t,y)
    $ x! E$ `& C9 h9 M. L. P5 w& ]) H
  6.                         dy = zeros(2,1);7 W6 K% k& n8 v; @& }0 L, C
  7.                         dy(1) = y(2)-f(t);8 m7 [4 U0 u  A) P: e3 O! s* B" P$ O
  8.                         dy(2) = y(1)*g(t)-y(2);1 I" ]0 C# X# F- @* V! k
  9.             end
    , v% A2 o& g! c" }
  10.             options = odeset('RelTol',1e-4,'AbsTol',[1e-5 1e-5]);
    : \+ z3 U2 d2 r3 `( Y/ X+ ~% J
  11.             %[T, y] = ode45(@rigid, [0 20], [1 2],options). _/ X6 G9 S8 L; _; ]" i
  12.             sol = ode45(@rigid, [0 20], [1 2],options);
    % |6 j# C; H7 L  P& n6 H- O
  13.             x=linspace(0,20,200000);
    ( k  P7 l" {+ r( i
  14.             y=deval(sol,x);, C; ~, l( ?3 I
  15.             res=y(1,:)+y(2,:);6 l2 o' D2 ~* l5 [( B5 I: E
  16.             idx=find(abs(res-0)<1e-4) %相加,当和小于误差运行范围的时候可以认为它就为0
    + g# q! ?, O) k+ J+ A, ^. g
  17.             xx=x(idx) %算出在x中的下标
    2 n  v* K, F8 |: j: I7 b6 ?
  18.             F=@(t)(f(t)+g(t));
    ) |+ T% n2 x. a
  19.             r=[]; %得出的解反代入方程求值,得到的值保存在r中。如果解正确,r中的值应该非常接近0: S8 Q8 o6 }# h( h
  20.             for i=1:size(xx,2)1 g; W. n2 g: {$ T/ A# L2 h" }; Z
  21.                     r = [r F(xx(i))]; %将解的值依次带入。当前问题在于r中的值都和0差的很远
    # M: n; |- g) U  {
  22.             end
    * g( C$ r! O) ?& j2 L4 \: b
  23.             r
    0 ?- Z) x$ T6 m  s% R0 ^9 ?
  24. end
复制代码
问题就是最后算出来的解再带入方程进行检验得出的结果和0差的好远,都到了1.73几。
$ w, y" G% [: n我也不知道哪里用错了,但猜测可能是由于微分方程里出现了f(t),g(t)的缘故
( z$ p: T  S4 t# s4 j
) B4 g' O  Y- \: O& U( Q求教大家我的程序是哪里出现了问题,该如何改呢?9 |! j  _: C- y) R0 h
谢谢( O0 w9 J& q3 R4 i: r% ]
- k5 B: s7 W' C6 k7 E  C
6 O* Y4 n+ P7 ]

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

问题

问题


作者: yzh07137    时间: 2015-8-4 20:45
算了,就知道这种论坛一般也没人0 @: R/ ?4 `/ C) J. ^





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