计算方法作业

1、 计算方法上机作业 机制1106班 U2011110763 陈成一、方程求根1.用牛顿迭代法求解下列方程的正根方程一：先画出f(x)= 的图像，大致确定根的大小，并由此确定初始点。M文件：function x=Newt_n(f_name,x0) % 创建Newton迭代的M文件x=x0; %迭代初值取x0xb=x+1; k=0;n=100;最大迭代次数del_x=0.01;求导数所用的间隔while abs(x-xb)0.000001 %打到要求误差，循环停止k=k+1;xb=x;if k=n break;end;y=feval(f_name,x);y_driv=(feval(f_name,x

2、+del_x)-y)/del_x;%用，间斜率代替导数x=xb-y/y_driv;%Newton迭代公式fprintf(k=%3.0f,x=%12.5e,y=%12.5e,yd=%12.5en,k,x,y,y_driv);end;fprintf(n Final answer=%12.6en,x);function y=eqn_1(x)y=log(x+1)-x.2;Newt_n(eqn_1,-0.5)k= 1, x=-1.82470e-001,y=-9.43147e-001,yd=2.97026e+000k= 2, x=-3.30076e-002,y=-2.34763e-001,yd=1.570

3、72e+000k= 3, x=-1.06356e-003,y=-3.46542e-002,yd=1.08484e+000k= 4, x=1.44016e-005,y=-1.06525e-003,yd=9.88214e-001k= 5, x=-2.19146e-007,y=1.44013e-005,yd=9.84990e-001k= 6, x=3.32970e-009,y=-2.19146e-007,yd=9.85034e-001 Final answer=3.329696e-009ans = 3.3297e-009 %也即是零Newt_n(eqn_1,4)k= 1, x=2.15747e+00

4、0,y=-1.43906e+001,yd=-7.81020e+000k= 2, x=1.28315e+000,y=-3.50489e+000,yd=-4.00872e+000k= 3, x=8.99412e-001,y=-8.20918e-001,yd=-2.13926e+000k= 4, x=7.69012e-001,y=-1.67397e-001,yd=-1.28373e+000k= 5, x=7.47720e-001,y=-2.09584e-002,yd=-9.84329e-001k= 6, x=7.46893e-001,y=-7.73188e-004,yd=-9.34897e-001k

5、= 7, x=7.46882e-001,y=-1.04149e-005,yd=-9.32973e-001k= 8, x=7.46882e-001,y=-1.29997e-007,yd=-9.32947e-001 Final answer=7.468817e-001ans =0.7469最终结果为：0.7469方程二：还是先画出图形，确定根的大致位置，然后根据图像中曲线走势找到收敛的初始点。而M文件可以沿用上面的。function y=eqn_1(x)y=exp(x)-5.*x.2; Newt_n(eqn_1,4)k= 1,x=5.71379e+000,y=-2.54018e+001,yd=1.

6、48221e+001k= 2,x=5.14868e+000,y=1.39780e+002,yd=2.47349e+002k= 3,x=4.82235e+000,y=3.96590e+001,yd=1.21530e+002k= 4,x=4.71816e+000,y=7.98116e+000,yd=7.66062e+001k= 5,x=4.70810e+000,y=6.57002e-001,yd=6.52923e+001k= 6,x=4.70794e+000,y=1.02918e-002,yd=6.42664e+001k= 7,x=4.70794e+000,y=8.23349e-005,yd=6.

7、42501e+001k= 8,x=4.70794e+000,y=6.48470e-007,yd=6.42500e+001 Final answer=4.707938e+000ans =4.7079Newt_n(eqn_1,1)k= 1,x=6.88208e-001,y=-2.28172e+000,yd=-7.31808e+000k= 2,x=6.11564e-001,y=-3.78006e-001,yd=-4.93195e+000k= 3,x=6.05364e-001,y=-2.67399e-002,yd=-4.31308e+000k= 4,x=6.05268e-001,y=-4.09485e

8、-004,yd=-4.26253e+000k= 5,x=6.05267e-001,y=-3.95813e-006,yd=-4.26175e+000 Final answer=6.052671e-001ans= 0.6053最终结果：4.7079 0.6053 方程三： 思路同上function y=eqn_1(x)y=x.3+2.*x-1; Newt_n(eqn_1,1)k= 1,x=6.02394e-001,y=2.00000e+000,yd=5.03010e+000k= 2,x=4.66118e-001,y=4.23383e-001,yd=3.10681e+000k= 3,x=4.5354

9、9e-001,y=3.35069e-002,yd=2.66588e+000k= 4,x=4.53398e-001,y=3.95932e-004,yd=2.63083e+000k= 5,x=4.53398e-001,y=2.09360e-006,yd=2.63041e+000 Final answer=4.533977e-001ans= 0.4534最终结果：0.4534 方程四：function y=eqn_1(x)y=(x+2).(0.5)-x;Newt_n(eqn_1,5)k= 1,x=2.09741e+000,y=-2.35425e+000,yd=-8.11085e-001k= 2,x=

10、2.00021e+000,y=-7.32032e-002,yd=-7.53140e-001k= 3,x=2.00000e+000,y=-1.58727e-004,yd=-7.50163e-001k= 4,x=2.00000e+000,y=-3.37165e-008,yd=-7.50156e-001 Final answer=2.000000e+000ans = 2.0000最终结果：2.00002.先用图解法确定初始点，然后再求方程 x sin x 1 = 0 x 0,10的所有根。必须先画出图像，确定根的个数，选择合适的初值点，否则极易得不到规定范围内的正确解。而且初值点选择很重要。因此图像

11、有求解的优点。Function y=eqn_1(x)y=x.*sin(x)-1;Newt_n(eqn_1,0.5)k= 1,x=1.32126e+000,y=-7.60287e-001,yd=9.25763e-001k= 2,x=1.10417e+000,y=2.80330e-001,yd=1.29134e+000k= 3,x=1.11416e+000,y=-1.38761e-002,yd=1.38935e+000k= 4,x=1.11416e+000,y=6.72563e-009,yd=1.38817e+000 Final answer=1.114157e+000ans =1.1142New

12、t_n(eqn_1,3)k= 1,x=2.79702e+000,y=-5.76640e-001,yd=-2.84083e+000k= 2,x=2.77312e+000,y=-5.51756e-002,yd=-2.30892e+000k= 3,x=2.77261e+000,y=-1.14804e-003,yd=-2.24109e+000k= 4,x=2.77260e+000,y=-7.70019e-006,yd=-2.23963e+000k= 5,x=2.77260e+000,y=-4.91882e-008,yd=-2.23962e+000Final answer=2.772605e+000an

13、s = 2.7726Newt_n(eqn_1,6)k= 1,x=6.48668e+000,y=-2.67649e+000,yd=5.49951e+000k= 2,x=6.43927e+000,y=3.10904e-001,yd=6.55805e+000k= 3,x=6.43912e+000,y=9.99083e-004,yd=6.52120e+000k= 4,x=6.43912e+000,y=7.40628e-007,yd=6.52106e+000 Final answer=6.439117e+000ans = 6.4391Newt_n(eqn_1,10)k= 1,x=9.27766e+000

14、,y=-6.44021e+000,yd=-8.91576e+000k= 2,x=9.31745e+000,y=3.59994e-001,yd=-9.04740e+000k= 3,x=9.31724e+000,y=-1.89206e-003,yd=-9.17150e+000k= 4,x=9.31724e+000,y=-3.11342e-006,yd=-9.17089e+000 Final answer=9.317243e+000ans = 9.3172最终结果：1.1142 2.7726 6.4391 6.4391 二、线性方程组1.计算下列矩阵的逆矩阵，并验证之。（1） 求 , B , BB=

15、1 4 5;2 1 2;8 1 1B = 1 4 5 2 1 2 8 1 1 inv(B)ans = -0.0400 0.0400 0.1200 0.5600 -1.5600 0.3200 -0.2400 1.2400 -0.2800 B*inv(B)ans = 1.0000 -0.0000 -0.0000 -0.0000 1.0000 -0.0000 -0.0000 -0.0000 1.0000 inv(B)*Bans = 1.0000 -0.0000 -0.0000 -0.0000 1.0000 -0.0000 -0.0000 0 1.0000(2) ， 求 ，C ，CC=0 5 1;-1

16、 6 3;8 -9 5C = 0 5 1 -1 6 3 8 -9 5 inv(C)ans = 0.5377 -0.3208 0.0849 0.2736 -0.0755 -0.0094 -0.3679 0.3774 0.0472 C*inv(C)ans = 1.0000 -0.0000 -0.0000 -0.0000 1.0000 0.0000 -0.0000 -0.0000 1.0000 inv(C)*Cans = 1.0000 0.0000 -0.0000 0 1.0000 0.0000 0 -0.0000 1.00002.用两种方法求解下列线性方程组（1）=（2）= 调用x=Ab 命令（1

17、）x= （2）x= (1)A=2 -1 0;-1 2 -1;0 -1 2:A = 2 -1 0 -1 2 -1 0 -1 2 b=1;2;3; x=Abx = 2.5000 4.0000 (方程一的解)3.5000(2)A=2 -1 1;-3 4 -1;0 -1 1A = 2 -1 1 -3 4 -1 0 -1 1 b=4;5;6; x=Abx = -1.0000 2.6667 (方程二的解) 8.6667 利用LU 分解(1) A=2 -1 0;-1 2 -1;0 -1 2A = 2 -1 0 -1 2 -1 0 -1 2 L U P=lu(A)L = 1.0000 0 0 -0.5000

18、1.0000 0 0 -0.6667 1.0000U = 2.0000 -1.0000 0 0 1.5000 -1.0000 0 0 1.3333P = 1 0 0 0 1 0 0 0 1 b=1;2;3; y=L(P*b)y = 1.0000 2.5000 4.6667 x=Uyx = 2.5000 4.0000 (方程一的解) 3.5000（2）A=2 -1 1;-3 4 -1;0 -1 1A = 2 -1 1 -3 4 -1 0 -1 1 L U P=lu(A)L = 1.0000 0 0 -0.6667 1.0000 0 0 -0.6000 1.0000U = -3.0000 4.00

19、00 -1.0000 0 1.6667 0.3333 0 0 1.2000P = 0 1 0 1 0 0 0 0 1 b=4;5;6； y=L(P*b)y = 5.0000 7.3333 10.4000 x=Uyx = -1.0000 2.6667 （方程二的解） 8.6667三、插值与拟合1.以下表格数据由函数f(x)= 得到 采用Lagrange 插值，求xi=0.2，0.6，1.0 处的函数值yi。以及误差值f(xi)-yi。xi=0.2 yi= f(xi)-yi=xi=0.6 yi= f(xi)-yi=xi=1.0 yi= f(xi)-yi=Lagrange插值法：function f

20、i=lagrange(x,y,xi)fi=zeros(size(xi);npl=length(y);for i=1:nplz=ones(size(xi);for j=1:nplif i=j,z=z.*(xi-x(j)/(x(i)-x(j);end;fi=fi+z*y(i);endx=0 0.4 0.8 1.2;y=1.0 1.491 2.225 3.320; xi=0.2 0.6 1.0;yi=lagrange(x,y,xi)yi =1.2225 1.8203 2.7200exp(xi)-yi = -0.0011 0.0019 -0.00172.已知飞机降落曲线为 ，式中y表示飞机高度，x 表

21、示飞机距指挥塔的距离。该函数满足条件 y (0) = 0、y(12000) =1000、y(0) = 0、y(12000) = 01) 试利用所满足的条件确定飞机降落曲线；2) 绘制出飞机降落曲线。列出关于四个未知数的方程，然后利用矩阵求解。由于未知数前系数过大，因此应该选择千米（km）作单位。A=0 0 0 1;12.3 12.2 12 1;0 0 1 0;3.*12.2 2.*12 1 0A = 0 0 0 1 1728 144 12 1 0 0 1 0 432 24 1 0 b=0;1;0;0； C=AbC = -0.0012 0.0208 0 0 x=0:0.001:12; plot(

22、x,-0.0012.*x.3+0.0208.*x.2)3.用三次多项式拟合下面数据，并做出图形x=0 0.2 0.4 0.6 0.8 1;y=0 7.78 10.68 8.37 3.97 0; a=polyfit(x,y,3)a = 41.5625 -101.6071 60.0768 -0.1179x=0:0.0001:1; plot(x, 41.5625.*x.3 -101.6071.*x.2+60.0768.*x-0.1179)四、数值积分1.用复合梯形求积法计算下列积分，取n=2、4、8、16（1） = = = =(2) = = = = function t=trapezia(a,b,n

23、,f_name)h=(b-a)/n;等分n等份fa=feval(f_name,a);左端点函数值fb=feval(f_name,b);右端点函数值t1=0;for k=1:n-1 x(k)=a+k*h; t1=t1+feval(f_name,x(k);end;t=h*(fa+2*t1+fb)/2;第k个小梯形的面积+前面的梯形面积function y=func_ts1(x)y=exp(x);endfunction y=func_ts2(x)y=1/(2+x);enda=2 4 8 16; for i=1:4;t1=trapezia(0,1,a(i),func_ts1)t2=trapezia(0

24、,1,a(i),func_ts2)endt 1= 1.7539t2 = 0.4083t 1= 1.7272t 2= 0.4062t 1= 1.7205t 2= 0.4056t1 = 1.7188t 2= 0.4055最终结果：2.用复合Simpson 求积法计算下列积分，取n=4、8、16、32 （1） = = = = （2） = = = = （3） = = = = M文件：function s=simpson(a,b,n,f_name)h=(b-a)/n;for k=1:2*n; x(k)=a+k*h/2;end;fa=feval(f_name,a);fb=feval(f_name,b);s

25、1=0;s2=0;for k=1:n s1=s1+ feval(f_name,x(2*k-1); s2=s2+ feval(f_name,x(2*k);end;s=h*(0.5*(fa-fb)+2*s1+s2)/3; function y=func_ts1(x)y=1/(2+cos(x);end /ts1function y=func_ts2(x)y=log(x+1)./x;end /ts2function y=func_ts3(x)y=1/(1+sin(x).2);end /ts3a=4 8 16 32;for i=1:4s1=simpson(0,pi,a(i),func_ts1)s2=si

26、mpson(10.-8,2,a(i),func_ts2)s3=simpson(0,pi/2,a(i),func_ts3)ends1 = 1.8138s2 = 1.4368s3 = 1.1107s1 = 1.8138s2 = 1.4367s3 = 1.1107s1 = 1.8138s2 = 1.4367s3 = 1.1107s1 = 1.8138s2 = 1.4367s3 = 1.1107最终结果：五、常微分方程1.试用Euler 法及改进Euler 法计算初值问题 取步长h=0.2。Euler法：function f=descent(t,y)f=y-2*t/y;a=0;=1;h=0.2;n=(

27、b-a)/h;t(1)=0; %t初值取零y(1)=1; %y对应初值取1for i=2:n+1t(i)=t(1)+(i-1)*h;y(i)=y(i-1)+h*descent(t(i-1),y(i-1);显示Euler法end;for i=1:n+1fprintf(Time=%f,t(i);fprintf( y=%fn,y(i);end;plot(t,y);Time=0.000000 y=1.000000Time=0.200000 y=1.200000Time=0.400000 y=1.373333Time=0.600000 y=1.531495Time=0.800000 y=1.681085Time=1.000000 y=1.826948改进Euler法：function f=algeb1(t,y)f=y-2*t/y;a=0;b=1;h=0.2;n=(b-a)/h;t(1)=0;y(1)