数学建模社区-数学中国
标题:
matlab求解微分方程结果不对
[打印本页]
作者:
yzh07137
时间:
2015-7-30 20:35
标题:
matlab求解微分方程结果不对
问题是这样的(见附图)
- B" u) ?% ?5 \4 m9 h
, H! m% k% `. e5 O7 N
然后我的程序如下
function Untitled
+ R v; W3 V3 |1 O& t' L
clear all;clc;
- p9 E0 C# D: l4 y
f=@(t)(2*sin(t)*(t<(4*pi)) + 0);
$ g- A. d9 J- [3 z9 R3 X( [8 _
g=@(t)(0+cos(t).*(t>=((7*pi)/2)));
; F' V0 g# f" z; V2 [1 a
function dy = rigid(t,y)
/ d# V' h) x' g' [' b: c( g3 C7 Q
dy = zeros(2,1);
! H$ W d% T0 \' j/ \9 [
dy(1) = y(2)-f(t);
% ]* q* t! `6 Q7 w w4 F
dy(2) = y(1)*g(t)-y(2);
/ x- j* M, i9 W8 h
end
8 h& x+ l/ |. R
options = odeset('RelTol',1e-4,'AbsTol',[1e-5 1e-5]);
5 D5 Y9 P' `' f6 k4 ~, E
%[T, y] = ode45(@rigid, [0 20], [1 2],options)
' B8 J0 v2 V- y+ M9 H W
sol = ode45(@rigid, [0 20], [1 2],options);
. k7 c- z$ N* o
x=linspace(0,20,200000);
: g/ C0 S; j- y# L- j+ C6 U3 s# a/ l
y=deval(sol,x);
# T( \& f2 v7 o
res=y(1,:)+y(2,:);
: s" n+ I, ^4 k, i
idx=find(abs(res-0)<1e-4) %相加,当和小于误差运行范围的时候可以认为它就为0
% a0 `( J+ T9 R( K
xx=x(idx) %算出在x中的下标
% k, l/ d- O# V6 z E' u3 I. }
F=@(t)(f(t)+g(t));
+ g7 ~. e# V; d+ A
r=[]; %得出的解反代入方程求值,得到的值保存在r中。如果解正确,r中的值应该非常接近0
6 q. R( q% @/ Q& L9 b7 ^0 A
for i=1:size(xx,2)
$ K" R: |4 a( N4 f8 G5 d- F3 b
r = [r F(xx(i))]; %将解的值依次带入。当前问题在于r中的值都和0差的很远
* S" J) ~/ k$ @3 Q5 S
end
! D [: Z! C0 D% L
r
/ h7 e* _4 X4 H
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)
2015-7-30 20:35 上传
点击文件名下载附件
问题
作者:
yzh07137
时间:
2015-8-4 20:45
算了,就知道这种论坛一般也没人
1 ]4 u, H- p7 ~
欢迎光临 数学建模社区-数学中国 (http://www.madio.net/)
Powered by Discuz! X2.5