常微分方程组(边值)

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

打靶法Shooting Method （shooting.m ）

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

function [x,a,b,n]=shooting(fun,x0,xn,eps)

if nargin<3

eps=1e-3;

end

x1=x0+rand;

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

c0=b(length(b),1);

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

c1=b(length(b),1);

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

n=1;

while (norm(c1-xn)>=eps & norm(x2-x1)>=eps)

x0=x1;x1=x2;

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

c0=b(length(b),1);

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

c1=b(length(b),1)

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

n=n+1;

end

x=x2;

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

()()⎪⎪⎩

⎪⎪⎨⎧==-=010004822y y y dx

y d

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

()⎪⎪⎪⎪⎩

⎪⎪⎪⎪⎨⎧==-==t dx dy y y y dx dy y dx dy x 001122

1048

命令：

x0=[0:0.1:10];

y0=32*((cos(5)-1)/sin(5)*sin(x0/2)-cos(x0/2)+1); 真实解

plot(x0,y0,'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))

函数：（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]=shooting('odebvp',10,0,1e-3)

计算结果：（eps=0.001）

t=11.9524

plot(x,y(:,1))

x0=[0:1:10];

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

plot(x0,y0,’o’)

有限差分法Finite Difference Methods FDM （difference.m ）

同上例：

4822y dx y d -=⇒4822

11i i i i y h y y y -=+--+ 2121

842h y y h y i i i =+⎪⎪⎭

⎫ ⎝⎛--+- 若划分为10个区间，则： ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡--=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣

⎡⎪⎪⎭⎫ ⎝⎛--⎪⎪⎭⎫ ⎝⎛--⎪⎪⎭⎫ ⎝⎛--⎪⎪⎭⎫ ⎝⎛----08880842114211421142222212212222h h h h y y y y h h h h n n M M O O O

函数：（difference.m ）

%有限差分法求常微分方程的边值问题

function [x,y]=difference(x0,xn,y0,yn,n)

h=(xn-x0)/n;

a=eye(n-1)*(-(2-h^2/4));

for i=1:n-2

a(i,i+1)=1;

a(i+1,i)=1;

end

b=ones(n-1,1)*8*h^2;

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)=yy(i-1);

end

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

命令：

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

计算结果：

x0=[0:0.1:10];

y0=32*((cos(5)-1)/sin(5)*sin(x0/2)-cos(x0/2)+1); 真实解plot(x0,y0,'r')

hold on

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

plot(x,y,’.’)

[x,y]=difference(0,10,0,0,10);

plot(x,y,’--’)

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

plot(x,y,’-.’)