最新上海电力学院数值计算方法上机实习题

2017数值计算方法上机实习报告

学院：

专业：

班级：

姓名：

学号：

数值计算方法上机实习题

1． 设?+=1

05dx x

x I n

n ， （1） 由递推公式n

I I n n 1

51+

-=-，从0=0.1823I ,1824.00=I 出发，计算20I ； （2） 20=0I ，20=10000I , 用n

I I n n 51

5111+-=--，计算0I ；

（3） 分析结果的可靠性及产生此现象的原因(重点分析原因)。

解：（1）分别令I 0的近似值为0.1823、0.1824，MATLAB 程序如下：

I=0.1823; %题中的已知数据 for n=1:20;

I=(-5)*I+1/n; %由递推公式所得 end

fprintf('I20=%f\n',I)

M=0.1824; %与I 的计算结果形成对比 for i=1:20;

M=(-5)*M+1/i; %由递推公式所得 end

fprintf('M20=%f\n',M) %% 输出结果

I20=-2055816073.851284 M20=7480927090.212283

（2）分别令I 20的近似值为0、10000，MATLAB 程序如下：

I=0; %赋予I20的初始值 for n=0:19;

I=(-1/5)*I+1/(5*(20-n)); %由递推公式所得 end fprintf('I0=%f\n',I) M=10000; for i=0:19;

M=(-1/5)*M+1/(5*(20-i));%由递推公式所得 end

fprintf('M0=%f\n',M) %% 输出结果

I0=0.182322 M0=0.182322 （3）分析：

由输出结果可看出第一种算法为不稳定算法，第二种算法为稳定算法。 由于误差

*

000I I e -=

02211*1*11*555)(5)15(15e e e I I n

I n I I I e n

n n n n n n n n n ===-=+--+

-=-=------

第一种算法为正向迭代算法，每计算一步误差增长5倍，虽然所给的初始值很接近，随着n 的增大，误差也越来越大。

*

n n n I I e -=

n n

n n n n n n n e I I n I n I e I I e 51)(51)5151(5151***111=-=+--+-==-=--- n n e e )5

1

(0=

第二种算法为倒向迭代算法，每计算一步误差缩小5倍，虽然所给的初始值之间差很多，随着n 的增大，误差也越来越小。算法趋近稳定，收敛，可以选择这种算法。 2． 求方程0210=-+x e x

的近似根，要求41105-+?<-k k x x ，并比较计算量。

（1） 在[0，1]上用二分法； （2） 取初值00=x ，并用迭代10

21x k e x -=+；

（3） 加速迭代的结果；

（4） 取初值00=x ，并用牛顿迭代法； （5） 分析绝对误差。 解：

（1）利用二分法，MATLAB 程序如下：

%% 二分法程序 clear all clc a=0;b=1;

f=@(x)(exp(x)+10*x-2); %@是定义函数句柄的运算符 c=(a+b)/2;%取区间中点 i=0;%分割次数

while abs(f(c))>5*10^(-4) %判断f(x)的精度是否满足要求 if f(a)*f(c)<0 b=c;c=(a+b)/2; elseif f(b)*f(c)<0 a=c;c=(b+a)/2; end i=i+1; end

fprintf('二分法运算次数为%i\n',i) fprintf('二分法计算结果为%f\n',c) %% 输出结果

二分法运算次数为13

二分法计算结果为0.090515

（2）用题目中给出的迭代法，取初始值x(1)=0，并用迭代x(i)=(2-exp(x(i-1)))/10，MATLAB 程序如下：

%% 不动点迭代

clear all

clc

x0=0;

x=x0;

for k=1:10000 %规定迭代次数上限

y=(2-exp(x))/10; %迭代结果存到y中

if abs(x-y)<5*10^(-4)

fprintf('初始值为x0%i\n迭代次数为%i\n',x0,k);

break

end

x=y;

if k==10000;

fprintf('迭代次数超出上限%i\n',k);

end

end

fprintf('迭代法计算结果为%f\n',y);

%% 输出结果

初始值x0为0

迭代次数为4

迭代法计算结果为0.090513

（3）利用加速迭代法，MATLAB程序如下：

%% 加速迭代法

x=0;

f=@(x)(( 2-exp(x))./10);

y=f(x); y0=f(y);

y1=x-((y-x)^2)/(y0-2*y+x);

i=1;

while abs(y1-x)> 5*10^(-4)

x=y1;

y=f(x);

y0=f(y);

y1=x-((y-x)^2)/(y0-2*y+x);

i=i+1;

end

%% 输出结果

y1= 0.090525

i=2

（4）牛顿迭代法：

取初始值x=0，MA TLAB程序如下：

%% 清空环境变量

clear all clc

%% 牛顿迭代法 x=0;

y=x-(exp(x)+10*x-2)./(exp(x)+10); i=1;

while abs(x-y)>0.0005 x=y;

y=x-(exp(x)+10*x-2)./(exp(x)+10); i=i+1; end %% 输出结果

y= 0.090525 i=2

（5）分析绝对误差

根据方程解求得对应的绝对误差如下表所示：

通过上述表格可以看出，二分法运算了11次，不动点迭代方法运算了4次， Atiken 加速迭代法运算了2次，牛顿迭代法运算2次.比较绝对误差可以发现 Atiken 加速迭代和牛顿迭代的计算结果的绝对误差较小。下面就其原因进行分析： 我们知道1212()()g x g x L x x -≤-,其中L 为压缩常数，并且01L ≤<。

进行误差估计：10*1k

k L x x x x L

-≤

--，当L 较小时，收敛较快，反之，当L 很靠近1时，收敛很慢。若

1L ≥时，则迭代不收敛。

由()1g x L '≤< ，()1()

lim =(0,1)!p k k k

e g x C C p e p *+→∞=≠≥常数，收敛阶为p 。

二分法和不动点迭代为线性收敛；Atiken 迭代和牛顿迭代是平方收敛，是超线性收敛的。其

中，二分法L=0.5；用迭代10

21x

k e x -=+，

2()==0.27101010

x x e e e g x '??-'≤≈ ???，比0.5小，

因此收敛比二分法快。不动点迭代为线性收敛，Atiken 迭代速度与牛顿迭代速度都是超线性收敛的，因而收敛速度较快。相比之下，只有二分法的收敛速度较慢,绝对误差最大。 3．钢水包使用次数多以后，钢包的容积增大，数据如下：

试从中找出使用次数和容积之间的关系，计算均方差。（用ax b

y c x

+=

+拟合） 解：拟合曲线的模型是ax b y c x +=+，将原模型变为)()(b ax x c y +=+，采用非线性最小二程

法。按照最小二乘原理，应选取参数a ，b ，c 使得表达式∑=+-+=

15

1k 2

))()((Z b ax x c y k

k k 达到极小值。具体的方法就是对Z 关于a 、b 、c 求偏导数，并置这些偏导数等于0，相对应的方程组如下所示：

???

?

??

???=+-+-=??=+-+-=??=+-+-=??∑∑∑---0))()((20))()((20))()((2151151151k k k k k

k k k k k k k k k b ax y x c y c Z b ax y x c b Z b ax y x c x a Z 通过对上述的偏导数方程组进行整理，可以将之写成的b AX =的形式，利用 MA TLA 充当计算器的作用，可以进行X 的求解，对应的就是a 、b 、c 三个参数的值。下面的程序就是进行求解方程组、画拟合曲线以及求取均方差的程序，具体如下： %求拟合方程，画拟合曲线

x=[2 3 4 5 6 7 8 9 10 11 12 13 14 15 16];

y=[6.42 8.2 9.58 9.5 9.7 10 9.93 9.99 10.49 10.59 10.60 10.8 10.6 10.9 10.76];

A=[sum(x.^2),sum(x),-sum(x.*y); sum(x),15,-sum(y); sum(x.*y),sum(y),-sum(y.^2)]; B=[sum(x.^2.*y);sum(x.*y);sum(y.^2.*x)]; f=A\B;

a=f(1);b=f(2);c=f(3); z=linspace(2,50,48);

Y=(f(1)*z+f(2))./(f(3)+z); plot(x,y,'r*',z,Y ,'b-');

fprintf('a=%f\nb=%f\nc=%f\n',a,b,c);

fprintf('拟合出来的方程式为:\n(%7.4f+x)y=%7.4f+%7.4fx\n',f(3),f(2),f(1)); %求均方差 for i=1:15

y1(i)=(f(2)+f(1)*x(i))/(f(3)+x(i)); end c=0;

for i=1:15

c=c+(y(i)-y1(i))^2; end

jfc=sqrt(c/15);

fprintf('均方差为%f\n',jfc) %%运行结果： a=11.340048 b=-12.495325 c=-0.340291

拟合出来的方程式为:

(-0.3403+x)y=-12.4953+11.3400x 均方差为0.220812

对应的拟合曲线如下图：

分析总结：此方法利用线性方法求解非线性问题，避免了非线性误差较大的问题。

4．设??????????

??----------------=410100141010014101101410010141001014A ，??????

???

?

??--=625250b ，b x =A 分析下列迭代法的收敛性，并求42

110-+≤-k

k x x 的近似解及相应的迭代次数。

（1） JACOBI 迭代；

（2） GAUSS-SEIDEL 迭代；

（3） SOR 迭代（取0.1:0.1:1.9ω=，找到迭代步数最少的*ω）。 解：

（1）JACOBI 迭代

%% JACOBI函数

function tx=jacobi(A,b,imax,x0,tol)

%初始值x0，次数imax，精度tol

del=10^(-10);

tx=[x0];

n=length(x0);

for i=1:n

dg=A(i,i);

if abs(dg)

disp('对角元太小');%防止出现溢出现象

return

end

end

for k=1:imax %jacobi循环

for i=1:n

sm=b(i);

for j=1:n

if j~=i

sm=sm-A(i,j)*x0(j);

end

end

x(i)=sm/A(i,i);

%x(1)到x(n)即完成一次迭代

end

tx=[tx;x];%矩阵中又加一行

if norm(x-x0)

return

else

x0=x;

end

end

%% 主函数程序

clear all

clc

A=[4 -1 0 -1 0 0;-1 4 -1 0 -1 0;0 -1 4 -1 0 -1;-1 0 -1 4 -1 0;0 -1 0 -1 4 -1;0 0 -1 0 -1 4]; b=[0 5 -2 5 -2 6];

x0=[0 0 0 0 0 0];

imax=100;%迭代次数

tol=10^(-4);%精度

tx=jacobi(A,b,imax,x0,tol);

for j=1:length(tx)

fprintf('%d%f%f%f%f%f%f\n',j,tx(j,1),tx(j,2),tx(j,3),tx(j,4),tx(j,5),tx(j,6));

end

%% 运行结果

10.0000000.0000000.0000000.0000000.0000000.000000 20.0000001.250000-0.5000001.250000-0.5000001.500000 30.6250001.0000000.5000001.0000000.5000001.250000 40.5000001.6562500.3125001.6562500.3125001.750000 50.8281251.5312500.7656251.5312500.7656251.656250 60.7656251.8398440.6796881.8398440.6796881.882813 70.9199221.7812500.8906251.7812500.8906251.839844 80.8906251.9252930.8505861.9252930.8505861.945313 90.9626461.8979490.9489751.8979490.9489751.925293 100.9489751.9651490.9302981.9651490.9302981.974487 110.9825741.9523930.9761961.9523930.9761961.965149 120.9761961.9837420.9674841.9837420.9674841.988098 130.9918711.9777910.9888951.9777910.9888951.983742 140.9888951.9924150.9848311.9924150.9848311.994448 150.9962081.9896390.9948201.9896390.9948201.992415 160.9948201.9964620.9929231.9964620.9929231.997410 170.9982311.9951670.9975831.9951670.9975831.996462 180.9975831.9983490.9966991.9983490.9966991.998792 190.9991751.9977450.9988731.9977450.9988731.998349 200.9988731.9992300.9984601.9992300.9984601.999436 210.9996151.9989480.9994741.9989480.9994741.999230 220.9994741.9996410.9992821.9996410.9992821.999737 230.9998201.9995090.9997551.9995090.9997551.999641 240.9997551.9998320.9996651.9998320.9996651.999877 250.9999161.9997710.9998861.9997710.9998861.999832 260.9998861.9999220.9998441.9999220.9998441.999943 270.9999611.9998930.9999471.9998930.9999471.999922 280.9999471.9999640.9999271.9999640.9999271.999973 290.9999821.9999500.9999751.9999500.9999751.999964 （2）GAUSS-SEIDEL迭代

%% GAUSS-SEIDEL函数

function tx=gseidel(A,b,imax,x0,tol)

%初始值x0，次数imax，精度tol

del=10^-10; tx=[x0];n=length(x0);

%tx是个二维矩阵，存储着每一步迭代的结果

for i=1:n

dg=A(i,i);

if abs(dg)

disp('对角元太小');

return

end

end

for k=1:imax %G-S循环

x=x0;

for i=1:n

sm=b(i);

for j=1:n

if j~=i

sm=sm-A(i,j)*x(j);

end

end

x(i)=sm/A(i,i);

end

tx=[tx;x];

if norm(x-x0)

return

else

x0=x;

end

end

%% 主函数程序

clear all

clc

A=[4 -1 0 -1 0 0;-1 4 -1 0 -1 0;0 -1 4 -1 0 -1;-1 0 -1 4 -1 0;0 -1 0 -1 4 -1;0 0 -1 0 -1 4];

b=[0 5 -2 5 -2 6];

x0=[0 0 0 0 0 0];

imax=100;%迭代次数

tol=10^(-4);%精度

tx=gseidel(A,b,imax,x0,tol);

for j=1:length(tx)

fprintf('%d %f %f %f %f %f %f\n',j,tx(j,1),tx(j,2),tx(j,3),tx(j,4),tx(j,5),tx(j,6)) end

%% 运行结果

1 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000

2 0.000000 1.250000 -0.187500 1.203125 0.113281 1.481445

3 0.613281 1.384766 0.51733

4 1.560974 0.606796 1.781033

4 0.73643

5 1.715141 0.764287 1.776880 0.818263 1.895638

5 0.873005 1.863889 0.884102 1.893843 0.913342 1.949361

6 0.939433 1.934219 0.944356 1.949283 0.958216 1.975643

7 0.970875 1.968362 0.973322 1.975603 0.979902 1.988306

8 0.985991 1.984804 0.987178 1.988268 0.990344 1.994381

9 0.993268 1.992698 0.993837 1.994362 0.995360 1.997299

10 0.996765 1.996490 0.997038 1.997291 0.997770 1.998702

11 0.998445 1.998313 0.998577 1.998698 0.998928 1.999376

12 0.999253 1.999189 0.999316 1.999374 0.999485 1.999700

13 0.999641 1.999610 0.999671 1.999699 0.999752 1.999856

14 0.999827 1.999813 0.999842 1.999855 0.999881 1.999931

15 0.999917 1.999910 0.999924 1.999931 0.999943 1.999967

16 0.999960 1.999957 0.999964 1.999967 0.999973 1.999984

（3）SOR迭代

%% SOR函数

function tx=sor(A,b,imax,x0,tol,w)

%初始值x0，次数imax，精度tol

del=10^-10;

tx=[x0];

n=length(x0);

for i=1:n

dg=A(i,i);

if abs(dg)

disp('对角元太小');

return

end

end

for k=1:imax %SOR循环

x=x0;

for i=1:n

sm=b(i);

for j=1:n

if j~=i

sm=sm-A(i,j)*x(j);

%先高斯迭代

end

end

x(i)=sm/A(i,i);

x(i)=w*x(i)+(1-w)*x0(i);

%比较高斯与超松弛迭代公式，补上省缺的项

end

tx=[tx;x];

if norm(x-x0)

return

else

x0=x;

end

end

%% 主函数程序

clear all

clc

A=[4 -1 0 -1 0 0;-1 4 -1 0 -1 0;0 -1 4 -1 0 -1;-1 0 -1 4 -1 0;0 -1 0 -1 4 -1;0 0 -1 0 -1 4]; b=[0 5 -2 5 -2 6];

x0=[0 0 0 0 0 0];

imax=100;迭代次数

tol=10^(-4);%精度

w=1.2;%给定不同的w,w=0.1:0.1:1.9，找出使SOR 迭代法收敛速度最快 tx=sor(A,b,imax,x0,tol,w); for j=1:size(tx,1)

fprintf('%d %f%f%f%f%f%f\n',j,tx(j,1),tx(j,2),tx(j,3),tx(j,4),tx(j,5),tx(j,6)) end %% 运行结果

1 0.0000000.0000000.0000000.0000000.0000000.000000

2 0.0000001.500000-0.1500001.4550000.2865001.840950

3 0.8865001.5069000.8708551.8221560.8937021.961177

4 0.8214171.9744120.9531531.9360500.9827511.988536

5 1.0088551.9885450.9833092.0052650.9981531.996732

6 0.9963721.9956411.0026291.9980940.9975092.000695

7 0.9988462.0005670.9992811.9990721.0005991.999825

8 1.0001231.9998870.9997792.0003360.9998951.999937

9 1.0000421.9999371.0001071.9999460.9999672.000035 10 0.9999572.0000220.9999791.9999821.0000181.999992 11 1.0000101.9999980.9999962.0000110.9999971.999999

W 值的范围在0~2之间，SOR 迭代才会收敛，令w=0.1:0.1:1.9，找出w*，使得SOR 迭代的收敛速度最快。W=0.1时需要迭代101次；w=1.9时需要迭代101次；w=1时需要迭代16次；w=0.9时需要迭代20次；w=1.1时需要迭代13次；w=1.3时需要迭代13次；w=1.2时需要迭代11次。由这几次试代可得出w=1.2时迭代次数最少，SOR 迭代法收敛速度最快。

总结：

由于A 为严格对角占有矩阵，根据定理可知雅克比迭代和GS 迭代都收敛，SOR 迭代收敛的必要条件是0

（1）JACOBI 迭代是根据A 分解成上三角、下三角和对角阵，将线性方程组的求解归结为对角方程组求解过程的重复，思路简单易于编程。迭代次数比较多，收敛速度慢；

（2）GAUSS-SEIDEL 迭代，是在JACOBI 收敛的基础上将计算出来的未知量马上投入下一个迭代方程中使用，使得迭代的速度加快，效果更好。但GS 迭代与JACOBI 迭代的收敛域并不能相互包容，所以不能相互代替。但如果两者都收敛时，一般说GS 迭代比JACOBI 迭代的收敛速度快。

（3）SOR 迭代为了加快收敛速度，在GS 的基础上加入一修正参数即松弛因子w ，使得收敛速度更快。但选着不同的w ，收敛速度不一样，为了达到最好的效果，要选着最合适的松弛因子。

5．用逆幂迭代法求????

? ??=111123136A 最接近于11的特征值和特征向量，准确到3

10-。

解：

%% 逆幂迭代法 A=[6 3 1;3 2 1;1 1 1]; v(:,1)=[1;0;1]; tol=1e-3; p=11;

L(1)=max(v(:,1));

u(:,1)=v(:,1)/L(1); k=1;

B=(A-p*eye(3)); while k>0

v(:,k+1)=B\u(:,k); L(k+1)=max(v(:,k+1)); u(:,k+1)=v(:,k+1)/L(k+1); err(k+1)=abs(1/L(k+1)-1/L(k)); if err(k+1)

迭代次数k =10 特征向量 =

4.4361 2.4363 1.0000 特征值 =7.8724

通过运行的结果可以看出，迭代运算了10次，最终求出的接近11的近似特征值为7.8724，对应的近似特征向量为[4.4361 2.4363 1.0000]’。 6．用经典R-K 方法求解初值问题 （1）??

?-+-='++-='x x y y y x y y y sin 2cos 22sin 22212211

，]10,0[∈x ，

??

?==3

)0(2

)0(21y y ； （2）???-+-='++-='x x y y y x y y y sin 999cos 999999998sin 22212

211

，]10,0[∈x ，

??

?==3)0(2

)0(2

1y y 。 和精确解???+=+=--x

e x y x

e x y x

x cos 2)(sin 2)(21比较，进行误差分析得到结论，图形显示精确解和数值解。 解：

（1）MATLAB 程序如下：

%% 用经典R-K 方法求解初值问题 clc;clear;

f=@(x,y1,y2)(-2*y1+y2+2*sin(x)) g=@(x,y1,y2) (y1-2*y2+2*cos(x)-2*sin(x)); h=0.1;

y1(1)=2;y2(1)=3;x(1)=0;

for i=1:100;

K1=f(x(i),y1(i),y2(i));

L1=g(x(i),y1(i),y2(i));

K2=f(x(i)+0.5*h,y1(i)+0.5*h*K1,y2(i)+0.5*h*L1); L2=g(x(i)+0.5*h,y1(i)+0.5*h*K1,y2(i)+0.5*h*L1); K3=f(x(i)+0.5*h,y1(i)+0.5*h*K2,y2(i)+0.5*h*L2); L3=g(x(i)+0.5*h,y1(i)+0.5*h*K2,y2(i)+0.5*h*L2); K4=f(x(i)+h,y1(i)+h*K3,y2(i)+h*L3);

L4=g(x(i)+h,y1(i)+h*K3,y2(i)+h*L3);

x(i+1)=x(i)+h;

y1(i+1)=y1(i)+h*(1/6)*(K1+2*K2+2*K3+K4);

y2(i+1)=y2(i)+h*(1/6)*(L1+2*L2+2*L3+L4);

end

x=0:0.1:10;

for i=1:101;

Y1(i)=2*exp(-x(i))+sin(x(i));

Y2(i)=2*exp(-x(i))+cos(x(i));

end

c=y1-Y1;

d=y2-Y2;

subplot(2,1,1)

plot(x,y1,'r*',x,Y1,'b-','LineWidth',2)

C=max(abs(y1-Y1))

C1=max(abs(y2-Y2))

legend('y1龙格库塔近似解'，'Y1精确解');

title('第一问曲线1精确解和近似解的对比'); ylabel('y1曲线Y1曲线')

subplot(2,1,2)

plot(x,y2,'r*',x,Y2,'b-','LineWidth',2)

legend('y2龙格库塔近似解'，'Y2精确解');

title('第一问曲线2精确解和近似解的对比'); ylabel('y2曲线Y2曲线')

%% 图形结果

误差结果：

C =

7.3449e-06

C1 =

7.4367e-06

（2）MATLAB程序如下：

%% 用经典R-K方法求解初值问题

clc;clear;

f=@(x,y1,y2)(-2*y1+y2+2*sin(x));

g=@(x,y1,y2)(998*y1-999*y2+999*cos(x)-999*sin(x)); h=0.001;

y1(1)=2;y2(1)=3;x(1)=0;

for i=1:10000;

K1=f(x(i),y1(i),y2(i));

L1=g(x(i),y1(i),y2(i));

K2=f(x(i)+0.5*h,y1(i)+0.5*h*K1,y2(i)+0.5*h*L1);

L2=g(x(i)+0.5*h,y1(i)+0.5*h*K1,y2(i)+0.5*h*L1);

K3=f(x(i)+0.5*h,y1(i)+0.5*h*K2,y2(i)+0.5*h*L2);

L3=g(x(i)+0.5*h,y1(i)+0.5*h*K2,y2(i)+0.5*h*L2);

K4=f(x(i)+h,y1(i)+h*K3,y2(i)+h*L3);

L4=g(x(i)+h,y1(i)+h*K3,y2(i)+h*L3);

x(i+1)=x(i)+h;

y1(i+1)=y1(i)+h*(1/6)*(K1+2*K2+2*K3+K4);

y2(i+1)=y2(i)+h*(1/6)*(L1+2*L2+2*L3+L4); end

x=0:0.001:10;

Y1=2*exp(-x)+sin(x);

Y2=2*exp(-x)+cos(x);%??è·?a

subplot(211)

plot(x,y1,'r*',x,Y1,'b-','LineWidth',2)

legend('y1龙格库塔近似解'，'Y1精确解'); title('第二问曲线1精确解和近似解的对比'); ylabel('y1曲线Y1曲线')

subplot(2,1,2)

plot(x,y2,'r*',x,Y2,'b-','LineWidth',2)

legend('y2龙格库塔近似解'，'Y2精确解'); title('第二问曲线2精确解和近似解的对比'); ylabel('y2曲线Y2曲线')

C=max(abs(y1-Y1))

C1=max(abs(y2-Y2))

%% 图形结果

误差结果：

C =

2.4263e-11

C1 =

2.3547e-08

四阶龙格库塔法所求出的解的精确度很高，很接近精确解。证明该方法是一种有效的算法。考虑到方程组的刚性问题，为了保证算法的稳定性，必须将步长限制在较小的范围内。第一组数据的刚性比为3大于1，为病态方程组，要求其步长限制在较小的范围内。第二组数据的刚性比为1000也大于1，为病态方程组，故步长也要限制在较小范围内。 7．用有限差分法求解边值问题（h=0.1），并图形显示。

?

?

?==-=+-''1)1()1(0

)1(2y y y x y . 解：%% 有限差分法求解边值问题

%% 追赶算法 function x=chase(a,b,c,f)

%求解线性方程组Ax=f ，其中A 使三对角阵 %a 是矩阵A 的小对角线元素a(1)=0 %b 是矩阵A 的对角线元素 %c 是矩阵A 的上对角线元素 %f 是方程组的右端向量 N=length(f);

x=zeros(1,N);y=zeros(1,N); d=zeros(1,N);u= zeros(1,N); %预处理 d(1)=b(1); for i=1:N-1 u(i)=c(i)/d(i);

d(i+1)=b(i+1)-a(i+1)*u(i); end %追的过程 y(1)=f(1)/d(1); for i=2:N

y(i)=(f(i)-a(i)*y(i-1))/d(i); end %赶的过程 x(N)=y(N); for i=N-1:-1:1 x(i)=y(i)-u(i)*x(i+1); end

%% 主程序 h=0.1; N=(1-(-1))/h; x(1)=-1; x(N-1)=1; y(1)=1; y(N-1)=1;

X(i)=x(1)+i*h;

q(i)=1+(X(i)^2);

B(i)=2/(h^2)+q(i);

end

for i=1:N-2

C(i)=-1/(h^2);

end

g(1)=0+1/(h^2);

g(N-1)=0+1/(h^2);

for i=2:N-2

g(i)=0;

end

y=chase([0,C],B,[C,0],g)

y1=[1,y,1];

X=-1:h:1;

plot(X,y1);

title('有限差分法');

xlabel('输入值x');

ylabel('y曲线');

%% 输出结果及图形

y =

Columns 1 through 10

0.9096 0.8357 0.7754 0.7267 0.6879 0.6577 0.6352 0.6195 0.6103 0.6073

Columns 11 through 19

0.6103 0.6195 0.6352 0.6577 0.6879 0.7267 0.7754 0.8357 0.9096

总结：步长分别取万分之一和十分之一时对比结果如上图，采用万分之一步长时可以近似为精确解，采用十分之一步长时可以认为是近似解，步长越小，运算越多。