實驗任務

（一）常微分方程的符號計算和數值解法基本操作

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('數值解','精確解')