数学建模 两点边值问题的两种数值解法及代码

常微分方程组两点边值问题的数值解法

3)1(1

)0(0

4===-''y y y y 可化为微分方程组3)1(1

)0(41221==='='y y y y y y

方法一：配置法

Matlab 程序：

function bvcollation

clc

solinit = bvpinit(linspace(0,1,20),[100 600]);% sol = bvp4c(@twoode,@twobc,solinit);

x = linspace(0,1,20);

y = deval(sol,x);

y'

plot(x,y(1,:),x,y(2,:));

end

%微分方程组

function dydx = twoode(x,y)

dydx = [ y(2)

4*y(1)];

end

%边值条件

function res = twobc(ya,yb)

res = [ ya(1)-1

yb(1)-3];

end

运行结果：

1.0000 -0.4203

0.9834 -0.2117

0.9777 -0.0055

0.9828 0.2007

0.9988 0.4091

1.0259 0.6220

1.0644 0.8419

1.1147 1.0710

1.1774 1.3121

1.2531 1.5677

1.3427 1.8407

1.4472

2.1341

1.5678

2.4512

1.7057

2.7954

1.8626 3.1707

2.0401

3.5811

2.2402 4.0313

2.4652 4.5261

2.7175 5.0712

3.0000 5.6724

方法二：打靶法

程序：积分用ode45,搜索用二分法

function shoot_first_try

clc

clear all

a=0;

b=1;

a2=-1;

b2=0;

s=0;

tol=1.0*10e-6;

h1=Fun(a2)

h2=Fun(b2)

s=bisect(@Fun,a2,b2,tol)

[t1,y1]=ode45(@f,[a,b],[1,s])

plot(t1,y1(:,1),t1,y1(:,2))

function xc=bisect(f,a,b,tol)

if sign(f(a))*sign(f(b))>=0

error('f(a)*f(b)<0 not satisfied!') end

fa=f(a);

fb=f(b);

k=0;

while (b-a)/2>tol

c=(a+b)/2;

fc=f(c);

if fc==0

break

end

if sign(f(a))*sign(f(b))<0 b=c;fb=fc;

else

a=c;fa=fc;

end

end

xc=(a+b)/2;

function z=Fun(s)

a=0;b=1;yb=3;

h=linspace(a,b,20);

y0=[1;s];

[t,y]=ode45(@f,h,y0);

z=y(end,1)-yb;

function ydot=f(t,y)

%ydot=[0;0];

ydot(1)=y(2);

ydot(2)=4*y(1);

ydot=[ydot(1);ydot(2)];

运行结果：

1.0000 -0.5000

0.9969 -0.4749

0.9940 -0.4499

0.9913 -0.4250

0.9887 -0.4001

0.9799 -0.3017

0.9736 -0.2041

0.9697 -0.1069

0.9683 -0.0100

0.9692 0.0868

0.9726 0.1839

0.9784 0.2814

0.9867 0.3797

0.9974 0.4788

1.0106 0.5792

1.0264 0.6811

1.0447 0.7846

1.0656 0.8901

1.0892 0.9978

1.1155 1.1080

1.1446 1.2210

1.1766 1.3370

1.2115 1.4564

1.2495 1.5794

1.2905 1.7064

1.3348 1.8377 1.3825 1.9735 1.4335

2.1143 1.4882 2.2603 1.5466 2.4120 1.6089 2.5698 1.6751 2.7339 1.7456 2.9049 1.8204

3.0832 1.8998 3.2692

1.9840 3.4633

2.0731

3.6661 2.1674 3.8781 2.2671

4.0998 2.3724 4.3317 2.4837 4.5745 2.5711 4.7637 2.6621 4.9596 2.7569

5.1625 2.8555 5.3726