常微分方程数值方法

常微分方程数值方法

1、欧拉方法：1,,1,0),,(1-=+=+n k y t hf y y k k k k .

function E=euler(f,a,b,ya,n)

% Input - f is the function entered as a string 'f'

% - a and b are the left and right end points

% - ya is the initial condition y(a)

% - n is the number of steps

% Output - E=[T' Y'] where T is the vector of abscissas and

% Y is the vector of ordinates

h=(b-a)/n;

T=zeros(1,n+1);

Y=zeros(1,n+1);

T=a:h:b;

Y(1)=ya;

for j=1:n

Y(j+1)=Y(j)+h*feval(f,T(j),Y(j));

end

E=[T' Y'];

【例】 用欧拉方法求解区间]3,0[内的初值问题：1)0(,2'=-=y y

t y 。

f=inline('(t-y)/2','t','y');a=0;b=3;ya=1;n=12;

%n=3,6,12,24,48,96...

E=euler(f,a,b,ya,n),plot(E(:,1),E(:,2),'r*'),hold on

符号解：y=dsolve('Dy=(t-y)/2','y(0)=1')

h=(3-0)/12;t=0:h:3;y=eval(y);[t' y']

用图比较数值解：(f 为ode 函数文件)

ode45('f',[0,3],1)

2、休恩(Huen)方法(即改进Euler 方法)：

1

,,1,0)],,(,(),([211-=+++=++n k y t hf y t f y t f h

y y k k k k k k k k

function H=heun(f,a,b,ya,n)

% Input - f is the function entered as a string 'f'

% - a and b are the left and right end points

% - ya is the initial condition y(a)

% - n is the number of steps

% Output - H=[T' Y'] where T is the vector of abscissas and

% Y is the vector of ordinates

h=(b-a)/n;

T=zeros(1,n+1);

Y=zeros(1,n+1);

T=a:h:b;

Y(1)=ya;

for j=1:n

k1=feval(f,T(j),Y(j));

k2=feval(f,T(j+1),Y(j)+h*k1);

Y(j+1)=Y(j)+(h/2)*(k1+k2);

end

H=[T' Y'];

【例】 用休恩方法求解区间]3,0[内的初值问题：1)0(,2'=-=y y

t y 。

f=inline('(t-y)/2','t','y');a=0;b=3;ya=1;n=12;

H=heun(f,a,b,ya,n),plot(H(:,1),H(:,2),'r*'),hold on

用图比较数值解：(f 为ode 函数文件)

ode45('f',[0,3],1)

3、四阶龙格－库塔方法：)22(643211k k k k h

y y k k ++++=+

function R=rk4(f,a,b,ya,n)

% Input - f is the function entered as a string 'f'

% - a and b are the left and right end points

% - ya is the initial condition y(a)

% - n is the number of steps

% Output - R=[T' Y'] where T is the vector of abscissas and

% Y is the vector of ordinates

h=(b-a)/n;

T=zeros(1,n+1);

Y=zeros(1,n+1);

T=a:h:b;

Y(1)=ya;

for j=1:n

k1=h*feval(f,T(j),Y(j));

k2=h*feval(f,T(j)+h/2,Y(j)+k1/2);

k3=h*feval(f,T(j)+h/2,Y(j)+k2/2);

k4=h*feval(f,T(j)+h,Y(j)+k3);

Y(j+1)=Y(j)+(k1+2*k2+2*k3+k4)/6;

end

R=[T' Y'];

【例】 用4阶龙格－库塔方法求解区间]3,0[内的初值问题：1)0(,2'=-=y y

t y 。

f=inline('(t-y)/2','t','y');a=0;b=3;ya=1;n=3;

R=rk4(f,a,b,ya,n)

plot(R(:,1),R(:,2),'r*'),hold on

用图比较数值解：(f 为ode 函数文件)

ode45('f',[0,3],1)

4、阿当姆斯－巴什弗斯－摩尔顿方法：（预测－校正公式）

⎪⎪⎩⎪⎪⎨⎧++-+=+-+-+=+--+---))

,(9195(24)

5559379(241121123p x f f f f h y y f

f f f h y p k k k k k k k k k k k

function A=adams(f,T,Y)

% Input - f is the function entered as a string 'f'

% - T is the vector of abscissas

% - Y is the vector of ordinates

% Remark:The first four coordinates of T and

% Y must have starting values obtained with RK4