MATLAB experiment＜6＞

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)

%发现所用节点很少