常微分方程组(边值)

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

打靶法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;

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

()()⎪⎪⎩

⎪⎪⎨⎧==-=01000482

2y y y

dx

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

()⎪⎪⎪⎪⎩

⎪

⎪⎪⎪⎨⎧==-==t dx dy

y y y dx

dy y dx dy x 0011221

048 命令： 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 -=⇒482211i i i i y h

y y y -=+--+ 2121

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

⎫

⎝⎛--+- 若划分为10个区间，则：

⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡--=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣

⎡⎪⎪⎭⎫ ⎝⎛--⎪⎪⎭⎫ ⎝⎛--⎪⎪⎭

⎫ ⎝⎛--⎪⎪⎭

⎫ ⎝⎛----08880842114211

421142222212212

2

22h 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,’-.’)