数学实验“微分方程组边值问题数值算法(打靶法,有限差分法)”实验报告(内含matlab程序)

西京学数学软件实验任务书
实验二十七实验报告
一、实验名称：微分方程组边值问题数值算法（打靶法，有限差分法）。

二、实验目的：进一步熟悉微分方程组边值问题数值算法（打靶法，有限差分法）。

三、实验要求：运用Matlab/C/C++/Java/Maple/Mathematica 等其中一种语言完成程序设计。

四、实验原理：
1．打靶法：
对于线性边值问题
⎩⎨⎧==∈=+'+''βα)(,)(],[)()()(b y a y b a x x f y x q y x p y (1)
假设L 是一个微分算子使：()()Ly y p x y q x y '''=++
则可得到两个微分方程：
)(1x f Ly =，α=)(1a y ，0)(1='a y
⇔)()()(111
x f y x q y x p y =+'+''，α=)(1a y ，0)(1='a y (2) 02=Ly ，0)(2=a y ，1)(2
='a y ⇔0)()(222
=+'+''y x q y x p y ，0)(2=a y ，1)(2='a y (3) 方程（2），（3）是两个二阶初值问题.假设1y 是问题（2）
的解，2y 是问题（3）的解，且2()0y b ≠，则线性边值问题（1）的解为:1122()
()()()()y b y x y x y x y b β-=+ 。

2．有限差分法：
基本思想是把连续的定解区域用有限个离散点构成的网格来代替，这些离散点称作网格的节点；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似， 积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组 ， 解此方程组就可以得到原问题在离散点上的近似解。

然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解。

五、实验内容：
%线性打靶法
function
[k,X,Y,wucha,P]=xxdb(dydx1,dydx2,a,b,alpha,beta,h) n=fix((b-a)/h); X=zeros(n+1,1); CT1=[alpha,0];
Y=zeros(n+1,length(CT1)); Y1=zeros(n+1,length(CT1)); Y2=zeros(n+1,length(CT1));
X=a:h:b;
Y1(1,:)= CT1;
CT2=[0,1];Y2(1,:)= CT2;
for k=1:n
k1=feval(dydx1,X(k),Y1(k,:))
x2=X(k)+h/2;y2=Y1(k,:)'+k1*h/2;
k2=feval(dydx1,x2,y2);
k3=feval(dydx1,x2,Y1(k,:)'+k2*h/2);
k4=feval(dydx1, X(k)+h,Y1(k,:)'+k3*h);
Y1(k+1,:)=Y1(k,:)+h*(k1'+2*k2'+2*k3'+k4')/6,k=k+1;
end
u=Y1(:,1)
for k=1:n
k1=feval(dydx2,X(k),Y2(k,:))
x2=X(k)+h/2;y2=Y2(k,:)'+k1*h/2;
k2=feval(dydx2,x2,y2);
k3=feval(dydx2,x2,Y2(k,:)'+k2*h/2);
k4=feval(dydx2, X(k)+h,Y2(k,:)'+k3*h);
Y2(k+1,:)=Y2(k,:)+h*(k1'+2*k2'+2*k3'+k4')/6,k=k+1;
end
v=Y2(:,1)
Y=u+(beta-u(n+1))*v/v(n+1)
for k=2:n+1
wucha(k)=norm(Y(k)-Y(k-1)); k=k+1;
end
X=X(1:n+1);Y=Y(1:n+1,:);k=1:n+1;wucha=wucha(1:k,:);
P=[k',X',Y,wucha'];
plot(X,Y(:,1),'ro',X,Y1(:,1),'g*',X,Y2(:,1),'mp')
xlabel('轴\it x'); ylabel('轴\it y')
legend('是边值问题的数值解y(x)的曲线','是初值问题1的数值解u(x)的曲线', '是初值问题2的数值解v(x)的曲线')
title('用线性打靶法求线性边值问题的数值解的图形')
%有限差分法
function
[k,A,B1,X,Y,y,wucha,p]=yxcf(q1,q2,q3,a,b,alpha,beta,h) n=fix((b-a)/h); X=zeros(n+1,1);
Y=zeros(n+1,1); A1=zeros(n,n);
A2=zeros(n,n); A3=zeros(n,n); A=zeros(n,n);B= zeros(n,1);
for k=1:n
X=a:h:b;
k1(k)=feval(q1,X(k)); A1(k+1,k)=1+h*k1(k)/2;
k2(k)=feval(q2,X(k));
A2(k,k)=-2-(h.^2)*k2(k);
A3(k,k+1)= 1-h*k1(k)/2; k3(k)=feval(q3,X(k));
end
for k=2:n
B(k,1)=(h.^2)*k3(k);
end
B(1,1)=(h.^2)*k3(1)-(1+h*k1(1)/2)*alpha;
B(n-1,1)=(h.^2)*k3(n-1)-(1+h*k1(n-1)/2)*beta;
A=A1(1:n-1,1:n-1)+A2(1:n-1,1:n-1)+A3(1:n-1,1:n-1);
B1=B(1:n-1,1);
Y=A\B1;Y1=Y'; y=[alpha;Y;beta];
for k=2:n+1
wucha(k)=norm(y(k)-y(k-1)); k=k+1;
end
X=X(1:n+1); y=y(1:n+1,1); k=1:n+1;
wucha=wucha(1:k,:); plot(X,y(:,1),'mp')
xlabel('轴\it x'); ylabel('轴\it y'),legend('是边值问题的数值解y(x)的曲线')
title('用有限差分法求线性边值问题的数值解的图形'),
p=[k',X',y,wucha'];。