【matlab】常微分方程的数值解法

实验任务

（一）常微分方程的符号计算和数值解法基本操作

1、课上例题：1、2、3、4

（二）专题实验（梯形格式）

编写梯形公式的程序，要求：

1. 程序要有通用性，例如：

          function [t,y]=tixing(f,a,b,y0,N,mytol)

          % 求解形式为y’=f(t,y)的常微分方程的梯形公式

          % a,b为自变量的取值范围

          % y0为初始值

          % N表示对区间[a,b]剖分N等份

  % mytol表示允许误差, 要求默认值为1e-6，用来作为校正终止的条件

1. 以例题4中的方程作为测试方程，画出数值解和精确解的图形，并做图例。在完成例4测试后，也可以再找其它方程作为测试方程。

（三）专题实验（改进Euler格式）

编写梯形公式的程序，要求：

（1）程序要有通用性，例如：

          function [t,y]=GJ_Euler(f,a,b,y0,N)

          % 求解形式为y’=f(t,y)的常微分方程的改进Euler公式

          % a,b为自变量的取值范围

          % y0为初始值

          % N表示对区间[a,b]剖分N等份

1. 以例4中的方程作为测试方程，画出数值解和精确解的图形，并做图例。在完成例4测试后，也可以再找其它方程作为测试方程。

（四）专题实验（变形Euler格式）

编写梯形公式的程序，要求：

（1）程序要有通用性，例如：

          function [t,y]=BX_Euler(f,a,b,y0,N)

          % 求解形式为y’=f(t,y)的常微分方程的变形Euler公式

          % a,b为自变量的取值范围

          % y0为初始值

          % N表示对区间[a,b]剖分N等份

1. 以例4中的方程作为测试方程，画出数值解和精确解的图形，并做图例。在完成例4测试后，也可以再找其它方程作为测试方程。

（五）专题实验（边值问题的有限差分方法）

编写课上的算例1、算例2，要求：

其中pfun, qfun, rfun，是方程中的系数和右端函数p(x),q(x),r(x),一般程序中不用单字母作为函数名。ceshi 是对编写的函数进行测试。myFD是有限差分方法的主程序，可以写为下列形式：

function [x,y]=myFD(alpha,beta,a,b,N)

%求常微分方程边值问题y’’+p(x)y’+q(x)y=r(x), a<x<b

                         y(a)=alpha, y(b)=beta

% N表示对区间[a,b]剖分N等份

%返回值x为[a,b]的剖分节点，

%返回值y为剖分节点上的数值解

（2）对课程上的算例1和算例2进行测试，每一个算例都要求画出数值解和精确解的图形，并做图例。测试的代码写在函数文件ceshi里，大家可以用其它方式编写。

（六）自选实验（可以不做）

1、编写Euler公式、向后Euler公式的程序，要求同专题实验（三）；

2、在二阶Runge-Kutta方法中选定适当的参数，编写程序，要求同专题实验（三）

实验程序和实验结果

1、实验程序：

>> y=dsolve('Dy=y^2','x')

 

y =

 

           0

 -1/(C2 + x)

2、实验程序：

>> y=dsolve('x*D2y-3*Dy=x^2','y(1)=0','y(5)=0','x')

 

y =

 

(31*x^4)/468 - x^3/3 + 125/468

3、实验程序：

>> y=dsolve('Df=f+g','Dg=-f+g','f(0)=1','g(0)=2','t')

y =

  包含以下字段的 struct:

    g: [1×1 sym]

    f: [1×1 sym]

>> y.f

 

ans =

 

exp(t)*cos(t) + 2*exp(t)*sin(t)

 

>> y.g

 

ans =

 

2*exp(t)*cos(t) - exp(t)*sin(t)

4、实验程序：

[x,y]=ode23(@funst,[0,1],pi/2);

yy=(x+pi/2)./cos(x);

plot(x,y,'*',x,yy,'-')

legend('数值解','解析解')

function yp=funst(x,y)

yp=y*tan(x)+sec(x);

end

（二）专题实验（梯形格式）

1、实验程序：

 function [x,y]=tixing(f,a,b,y0,N,mytol)

  % 求解形式为y’=f(t,y)的常微分方程的梯形公式

  % a,b为自变量的取值范围

  % y0为初始值

  % N表示对区间[a,b]剖分N等份

  % mytol表示允许误差, 要求默认值为1e-6，用来作为校正终止的条件

  %nargin表示参数个数

  if nargin<=5

      mytol=1e-6;

  end

  h=(b-a)/N;

  x=a:h:b;

  y(1)=y0;

  for i=1:length(x)-1

      y(i+1)=y(i)+h*f(x(i),y(i));

      y1=y(i+1)+1;

      while abs(y(i+1)-y1)>=mytol

          y1=y(i+1);

          y(i+1)=y(i)+h/2*(f(x(i),y(i))+f(x(i+1),y1));

      end

  end

 end

1. 实验程序：

先建立函数文件funst.m

function yp=funst(x,y)

yp=y*tan(x)+sec(x);

end

求解程序：

[x,y]=tixing(@funst,0,1,pi/2,10);

yy=(x+pi/2)./cos(x);

plot(x,y,'*',x,yy,'-')

legend('数值解','解析解')

（三）专题实验（改进Euler格式）

1、实验程序：

function [x,y]=GJ_Euler(f,a,b,y0,N)

  % 求解形式为y’=f(t,y)的常微分方程的梯形公式

  % a,b为自变量的取值范围

  % y0为初始值

  % N表示对区间[a,b]剖分N等份

  % mytol表示允许误差, 要求默认值为1e-6，用来作为校正终止的条件

  %nargin表示参数个数

  if nargin<=5

      mytol=1e-6;

  else

      error

  end

  h=(b-a)/N;

  x=a:h:b;

  y(1)=y0;

  for i=1:length(x)-1

      y(i+1)=y(i)+h*f(x(i),y(i));

      y(i+1)=y(i)+h/2*(f(x(i),y(i))+f(x(i+1),y(i+1)));

  end

 end

2、实验程序：

先建立函数文件funst.m

function yp=funst(x,y)

yp=y*tan(x)+sec(x);

end

求解程序:

[x,y]=GJ_Euler(@funst,0,1,pi/2,10);

yy=(x+pi/2)./cos(x);

plot(x,y,'*',x,yy,'-')

legend('数值解','解析解')

• 专题实验（变形Euler格式）

1. 实验程序：

function [x,y]=BX_Euler(f,a,b,y0,N)

  % 求解形式为y’=f(t,y)的常微分方程的梯形公式

  % a,b为自变量的取值范围

  % y0为初始值

  % N表示对区间[a,b]剖分N等份

  % mytol表示允许误差, 要求默认值为1e-6，用来作为校正终止的条件

  %nargin表示参数个数

  if nargin<=5

      mytol=1e-6;

  else

      error

  end

  h=(b-a)/N;

  x=a:h:b;

  y(1)=y0;

  for i=1:length(x)-1

      y(i+1)=y(i)+h*f(x(i)+h/2,y(i)+h/2*f(x(i),y(i)));

  end

 end

1. 实验程序：

先建立函数文件funst.m

function yp=funst(x,y)

yp=y*tan(x)+sec(x);

end

求解程序：

[x,y]=BX_Euler(@funst,0,1,pi/2,10);

yy=(x+pi/2)./cos(x);

plot(x,y,'*',x,yy,'-')

legend('数值解','解析解')

（五）专题实验（边值问题的有限差分方法）

1、实验程序：

function [x,y]=myFD(alpha,beta,a,b,N)

%求常微分方程边值问题y''+p(x)y’+q(x)y=r(x),a<x<b,y(a)=alpha,y(b)=beta

% N表示对区间[a,b]剖分N等份

%返回值x为[a,b]的剖分节点

%返回值y为剖分节点上的数值解

h=(b-a)/N;

x=a:h:b;

y=zeros(length(x),1);

y(1)=alpha; y(length(x))=beta;

f=zeros(N-1,1);

A=zeros(N-1,N-1);

for i=1:N-1

    if i==1

        f(i,1)=h^2*rfun(x(2))-(1-h/2*pfun(x(2)))*alpha;

        A(1,1)=-2+h^2*qfun(x(2));

        A(1,2)=1+h/2*pfun(x(2));

    elseif i==N-1

        f(i,1)=h^2*rfun(x(N))-(1+h/2*pfun(x(N)))*beta;

        A(N-1,N-2)=1-h/2*pfun(x(N));

        A(N-1,N-1)=-2+h^2*qfun(x(N));

    else

        f(i,1)=h^2*rfun(x(i+1));

        A(i,i-1)=1-h/2*pfun(x(i+1));

        A(i,i)=-2+h^2*qfun(x(i+1));

        A(i,i+1)=1+h/2*pfun(x(i+1));

    end

end

y(2:N,1)=A\f;

y

end

2、实验程序：

function px=pfun(x)

px=3;

end

function qx=qfun(x)

qx=2;

end

function rx=rfun(x)

rx=sin(x)+3*cos(x);

end

[x,y]=myFD(0,1,0,pi/2,10);

yy=sin(x);

plot(x,y,'*',x,yy,'-')

legend('数值解','精确解')

1. 实验程序：

function px=pfun(x)

px=2+x^2;

end

function qx=qfun(x)

qx=1+2*x^2;

end

function rx=rfun(x)

rx=-sin(x)+(2+x^2)*cos(x)+(1+2*x^2)*sin(x);

end

[x,y]=myFD(0,1,0,pi/2,10);

yy=sin(x);

plot(x,y,'*',x,yy,'-')

legend('数值解','精确解')