数学建模社区-数学中国

标题: matlab求解微分方程结果不对 [打印本页]
作者: yzh07137    时间: 2015-7-30 20:35
标题: matlab求解微分方程结果不对
问题是这样的（见附图）+ {/ I* w- Z" w3 u* t$ D4 M

( T5 j7 n; F1 w8 e1 E0 }2 T+ ~# q然后我的程序如下
  1. function Untitled+ r- [! Z; s- g
  2.         clear all;clc;, o! B! Z! v9 L+ g2 @
  3.         f=@(t)(2*sin(t)*(t<(4*pi)) + 0);6 {! a* W* i+ Q$ {1 v
  4.         g=@(t)(0+cos(t).*(t>=((7*pi)/2)));
    8 G7 F9 Q. E4 N+ p+ c3 e
  5.             function dy = rigid(t,y)
    * ~& p7 e7 Y1 m' S* {# _: R) H- A
  6.                         dy = zeros(2,1);: |3 U' D4 f6 r, c3 ^
  7.                         dy(1) = y(2)-f(t);
    + p2 _% W* Z. ~. K2 A3 j
  8.                         dy(2) = y(1)*g(t)-y(2);
    $ ^) U0 V) ~$ B9 C# B
  9.             end
    - x+ \! r+ X# ]% \: B
  10.             options = odeset('RelTol',1e-4,'AbsTol',[1e-5 1e-5]);
    ' Y3 Q  ~5 W" r6 z6 z
  11.             %[T, y] = ode45(@rigid, [0 20], [1 2],options)
    9 E6 h4 p% Y. O9 p/ d" J
  12.             sol = ode45(@rigid, [0 20], [1 2],options);, ^  q5 [; ~6 B8 N- m
  13.             x=linspace(0,20,200000);8 ~$ ~1 F5 V$ x' Q) Y8 C6 d
  14.             y=deval(sol,x);' C/ y' q/ w+ j5 `- M5 n6 ?, j
  15.             res=y(1,:)+y(2,:);
    $ D) {8 o0 M% [* g+ M) H2 `. y
  16.             idx=find(abs(res-0)<1e-4) %相加，当和小于误差运行范围的时候可以认为它就为0
    0 ?% N; g( g& h
  17.             xx=x(idx) %算出在x中的下标
    ' M4 m* i. j. [7 X- h( b
  18.             F=@(t)(f(t)+g(t));' p! B1 C7 h$ p0 F9 [5 R
  19.             r=[]; %得出的解反代入方程求值，得到的值保存在r中。如果解正确，r中的值应该非常接近0
    ( M! ~/ V* s# V
  20.             for i=1:size(xx,2)
    3 s+ o5 ~# ~- C/ I6 J0 n6 d0 h; C
  21.                     r = [r F(xx(i))]; %将解的值依次带入。当前问题在于r中的值都和0差的很远
    # W* a6 p' p) @. x: E
  22.             end
    1 }8 K- u3 B2 U3 l5 j/ ^
  23.             r! Q$ n2 V0 p% _; @4 ~
  24. end
复制代码
问题就是最后算出来的解再带入方程进行检验得出的结果和0差的好远，都到了1.73几。
0 g7 f& S4 d  d; ]' B2 I' h! n我也不知道哪里用错了，但猜测可能是由于微分方程里出现了f(t)，g(t)的缘故, u3 _$ ~2 V1 X$ }- j" {0 f0 O
/ |2 D+ |9 K- f% h$ R! F
求教大家我的程序是哪里出现了问题，该如何改呢？$ ~  }& M5 U6 R; P
谢谢
5 X+ O& Z* o. S; d$ D* O8 [8 {: u! y
( }) Y  w% N1 u% o1 j3 r4 V

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

问题

问题

作者: yzh07137    时间: 2015-8-4 20:45
算了，就知道这种论坛一般也没人
5 S) G5 L' F8 x* t& X& a



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