常微分方程组(边值)

常微分方程组边值问题解法

打靶法Shooti ng Method (shoot in g.m )

% 丁靶法求常微分方程的边值问题

function [x,a,b ,n]=shooti ng(fu n, xO,x n, eps) if nargin<3

eps=1e-3;

end

x1=x0+ra nd;

[a,b]=ode45(fu n, [0,10],[0,x0]');

c0=b(le ngth(b),1);

[a,b]=ode45(fu n, [0,10],[0,x1]');

c1=b(le ngth(b),1);

x2=x1-(c1-x n)*(x1-x0)/(c1-c0);

n=1;

while (no rm(c1-x n)>=eps & no rm(x2-x1)>=eps) x0=x1;x 仁x2;

[a,b]=ode45(fu n,[ 0,10],[0,x0]');

cO=b(le ngth(b),1);

[a,b]=ode45(fu n,[ 0,10],[0,x1]');

c1= b(le ngth(b),1)

x2=x1-(c1-x n)*(x1-x0)/(c1-c0);

n=n+1;

end

x=x2;

应用打靶法求解下列边值问题:

y 10 0 解：将其转化为常微分方程组的初值问题

命令:

xO=[O:O.1:1O];

y0=32*((cos(5)-1)/si n( 5)*si n(x0/2)-cos(x0/2)+1); plot(xO,yO,'r')

hold on [x,y]=ode45('odebvp',[0,10],[0,2]'); plot(x,y(:,1))

[x,y]=ode45('odebvp',[0,10],[0,5]'); plot(x,y(:,1))

[x,y]=ode45('odebvp',[0,10],[0,8]'); plot(x,y(:,1)) [x,y]=ode45('odebvp',[0,10],[0,10]'); plot(x,y(:,1))

dy i dx y 2 dy 2 dx y i 0

y 4 y o

dy

dx X0

真实解

30

'

12^4567^9 10

函数：(odebvp.m)

泌值常微分方程(组)函数

function f=odebvp(x,y)

f(1)=y(2)；

f(2)=8-y(1)/4;

f=[f(1)；f(2)]；

命令：

[t,x,y, n]=shoot in g('odebvp',10,0,1e-3)

计算结果：(eps=0.001 )

t=11.9524

plot(x,y(:,1))

x0=[0:1:10];

y0=32*((cos(5)-1)/si n( 5)*si n(x0/2)-cos(x0/2)+1); hold on

plot(xO,yO, ' o')

有限差分法 Finite Differenee Methods FDM

同上例:

Y i i

2

Y i y i i

h 2 2 ——

1

4

函数：(differe nce.m )

%有限差分法求常微分方程的边值问题 function [x,y]=differe nce(xO,x n,yO,yn,n) h=(x n-xO)/n;

a=eye( n-1)*(-(2-h A 2/4)); for i=1: n-2 a(i,i+1)=1; a(i+1,i)=1; end

b=o nes( n-1,1)*8*hA2; b(1)=b(1)-0; b(n-1)=b( n-1)-0; yy=a\b;

x(1)=x0;y(1)=y0; for i=2: n x(i)=x0+(i-1)*h;

Y i i 2

若划分为10个区间，则：

h- Y i Y i 1 8h 2

4

(differe nce.m )

1

h 2 2 4

1

1

Y 1 Y 2

8h 2

8h 2

0 h 2

Y n 2 8h 2

2

4 1

Y n 1

8h 2 0

.2

h 1

2

4

y(i)=yy(i-i);

end

x(n)=xn;y(n)=yn;

命令：

[x,y]=differe nce(0,10,0,0,100);

计算结果：

xO=[O:O.1:1O];

y0=32*((cos(5)-1)/si n( 5)*si n(x0/2)-cos(x0/2)+1);

plot(xO,yO,'r')

hold on

[x,y]=differe nce(0,10,0,0,5);

plot(x,y,'.')

[x,y]=differe nce(0,10,0,0,10);

plot(x,y,'--')

[x,y]=differe nce(0,10,0,0,50);

plot(x,y,'-.')

正交配置法Orthogonal Collocatioin Methods CM

构造正交矩阵函数(collmatrix.m )

%E交配置矩阵(均用矩阵法求对称性与非对称性正交配置矩阵) function [am,bm,wm,an,bn,wn]=collmatrix(a,m,fm,n,fn) 真实解

的

3 4 6 7 9 9 1D

X