1,解非刚性微分方程
[t,y] = ode45(odefun,tspan,y0)
(其中 tspan = [t0 tf]
)求微分方程组 y′=f(t,y) 从 t0
到 tf
的积分,初始条件为 y0
。解数组 y
中的每一行都与列向量 t
中返回的值相对应。
1. 解下列微分方程。
%Excercise 1(1)
fun=inline('x+y','x','y');
[t,y]=ode45(fun,[0 1 2 3],1) %注意由于初值为y(0)=1, [0 1 2 3]中0不可缺
%Excercise 1(2)
fun=inline('[-2*y(1)-3*y(2);2*y(1)+y(2)]','t','y');%原题有错,应改为-x'=2x+3y
[t,y]=ode45(fun,[0 10],[-2.7;2.8])
plot(y(:,1),y(:,2))
%Excercise 1(3)
%高阶导数y''化为一阶(y')',多变量(y,y')化为单变量x.
%令x(1)=y,x(2)=y',化为方程组
%x(1)'=x(2),x(2)'=0.01*x(2)^2-2*x(1)+sin(t)
%初始值x(1)=0,x(2)=1.
%运行下列指令
clear;close;
fun=@(t,x)[x(2);0.01*x(2)^2-2*x(1)+sin(t)];%fun表示两个方程的右端,注意第一个x(2)表示x(1)的导函数。
[t,x]=ode45(fun,[0 5],[0;1]);x(end,1)
plot(t,x(:,1))
3. 求解刚性方程组
dfun=inline('[-1000.25*y(1)+999.75*y(2)+0.5;999.75*y(1)-1000.25*y(2)+0.5]','x','y');
[x,y]=ode45(dfun,[0,50],[1;-1]);length(x)
%发现所用节点很多
[x,y]=ode15s(dfun,[0,50],[1;-1]);length(x)
%发现所用节点很少