数学建模社区-数学中国

标题: matlab求解微分方程结果不对 [打印本页]
作者: yzh07137    时间: 2015-7-30 20:35
标题: matlab求解微分方程结果不对
问题是这样的（见附图）- B" u) ?% ?5 \4 m9 h

, H! m% k% `. e5 O7 N然后我的程序如下
  1. function Untitled
    + R  v; W3 V3 |1 O& t' L
  2.         clear all;clc;
    - p9 E0 C# D: l4 y
  3.         f=@(t)(2*sin(t)*(t<(4*pi)) + 0);
    $ g- A. d9 J- [3 z9 R3 X( [8 _
  4.         g=@(t)(0+cos(t).*(t>=((7*pi)/2)));; F' V0 g# f" z; V2 [1 a
  5.             function dy = rigid(t,y)
    / d# V' h) x' g' [' b: c( g3 C7 Q
  6.                         dy = zeros(2,1);
    ! H$ W  d% T0 \' j/ \9 [
  7.                         dy(1) = y(2)-f(t);
    % ]* q* t! `6 Q7 w  w4 F
  8.                         dy(2) = y(1)*g(t)-y(2);/ x- j* M, i9 W8 h
  9.             end
    8 h& x+ l/ |. R
  10.             options = odeset('RelTol',1e-4,'AbsTol',[1e-5 1e-5]);5 D5 Y9 P' `' f6 k4 ~, E
  11.             %[T, y] = ode45(@rigid, [0 20], [1 2],options)
    ' B8 J0 v2 V- y+ M9 H  W
  12.             sol = ode45(@rigid, [0 20], [1 2],options);. k7 c- z$ N* o
  13.             x=linspace(0,20,200000);: g/ C0 S; j- y# L- j+ C6 U3 s# a/ l
  14.             y=deval(sol,x);# T( \& f2 v7 o
  15.             res=y(1,:)+y(2,:);
    : s" n+ I, ^4 k, i
  16.             idx=find(abs(res-0)<1e-4) %相加，当和小于误差运行范围的时候可以认为它就为0
    % a0 `( J+ T9 R( K
  17.             xx=x(idx) %算出在x中的下标% k, l/ d- O# V6 z  E' u3 I. }
  18.             F=@(t)(f(t)+g(t));
    + g7 ~. e# V; d+ A
  19.             r=[]; %得出的解反代入方程求值，得到的值保存在r中。如果解正确，r中的值应该非常接近0
    6 q. R( q% @/ Q& L9 b7 ^0 A
  20.             for i=1:size(xx,2)$ K" R: |4 a( N4 f8 G5 d- F3 b
  21.                     r = [r F(xx(i))]; %将解的值依次带入。当前问题在于r中的值都和0差的很远
    * S" J) ~/ k$ @3 Q5 S
  22.             end! D  [: Z! C0 D% L
  23.             r/ h7 e* _4 X4 H
  24. end
复制代码
问题就是最后算出来的解再带入方程进行检验得出的结果和0差的好远，都到了1.73几。
3 f7 v$ c2 O/ n7 |我也不知道哪里用错了，但猜测可能是由于微分方程里出现了f(t)，g(t)的缘故
. S- T$ Z1 q% M) f' K
7 f4 W, T% a# `( R" R( M! C求教大家我的程序是哪里出现了问题，该如何改呢？1 ?2 g+ e  c' x2 R( F% ^: _
谢谢
+ P2 l) M+ a! o( @* f# r) ?1 A/ B1 G8 }. i8 m8 b2 K5 B) ?0 M* ?

* \5 I, d2 M" E( j

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

问题

问题

作者: yzh07137    时间: 2015-8-4 20:45
算了，就知道这种论坛一般也没人1 ]4 u, H- p7 ~




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