MATLAB数学建模(8)-常微分方程求解

使用改进欧拉方法求解一般的非刚性常微分方程。

function [x,y] = eulerpro(fun,x0,xfinal,y0,n);
if nargin < 5,n = 50;end;
h = (xfinal - x0)/n;
x(1) = x0;y(1) = y0;
for i = 1:n
    x(i+1) = x(i) + h;
    y1 = y(i) + h*feval(fun,x(i),y(i));
    y2 = y(i) + h*feval(fun,x(i+1),y1);
    y(i+1) = (y1 + y2)/2;
end

function f = doty(x,y);
f = -2*y + 2*x^2 +2*x;

运行结果:

>> [x,y] = eulerpro('doty',0,0.5,1,10)

x =

  Columns 1 through 7

         0    0.0500    0.1000    0.1500    0.2000    0.2500    0.3000

  Columns 8 through 11

    0.3500    0.4000    0.4500    0.5000


y =

  Columns 1 through 7

    1.0000    0.9076    0.8293    0.7641    0.7112    0.6701    0.6400

  Columns 8 through 11

    0.6204    0.6107    0.6105    0.6194

随着迭代次数的增加，不断接近精确解。
在MATLAB的工具箱中，提供了几个解非刚性常微分方程的函数，如ode45,ode23,ode113,其中ode45采用四五阶RK方法，是解决这类问题的首选方法，ode23采用二三阶RK方法，ode113采用多步法，效率比ode45高。
如求解一阶常微分方程:
t^2 y’ = y + 3t,y(1) = -2,1

odefun=@(t,y) (y+3*t)/t^2; %定义函数
tspan=[1 4]; %求解区间
y0=-2; %初值
[t,y]=ode45(odefun,tspan,y0);
plot(t,y) %作图
title('t^2y''=y+3t,y(1)=-2,1<t<4')
legend('t^2y''=y+3t')
xlabel('t')
ylabel('y') % 精确解

结果: