您的当前位置:首页数值分析课程设计大作业

数值分析课程设计大作业

2020-12-20 来源:飒榕旅游知识分享网
课程设计

设计题目 《数值分析》课程设计

学生姓名 ****

学 号 ####

专业班级

指导教师

2013年07月20日

设计 《数值分析》课程设计 题目 实验一:三题全做;实验二:四题全做;实验三:3.1 3.2已成绩 课程设计主要内容 做;实验四 4.1 4.2 4.3 4.4已做;实验五:5.1 5.2 5.3 5.4已做;实验六:6.2 6.3 6.4 已做;实验七:7.3 7.5 7.6 7.8 已做;实验八:8.2 8.8已做。共计26道题。 部分题通过matlab中notebook实现程序,notebook使用非常方便。由于篇幅有限部分题调用的函数源代码没有给出,本作业所使用的m文件均在附录中 建议:从学生的工作态度、工作量、设计(论文)的创造性、学术性、实用性及书面表达能力等方面给出评价。 指导教师评语 签名: 20 年 月 日 1.1水手、猴子和椰子问题

算法分析:设椰子起初的数目为p0,第一至第五次猴子在夜里藏椰子后,椰子的数目分别

为p0,p1,p2,p3,p4,再设最后每个人分得x个椰子,由题意得:

41pk1(pk1),k0,1,2,3,4.x(p51)所以 p5=5x1

55利用逆向递推方法求解:

n=input('n='); for x=1:n p=5*x+1; for k=1:5 p=5*p/4+1; end

if p==fix(p) break end end

disp([x,p])

执行代码后得: n= 1023 15621 (输入n=1000000) 即最后每个人分得1023个椰子,椰子总数为15621

1.2当n0,1,2,由

1enxdx ,100时,选择稳定的算法计算积分Inx0e101ex10I110I0xdx1,0e10 (n1)xnx1e110e1nxnIn110Indxedx(1e)x00e10n 得

1I(1I1)010I1[1(1en)I],n100,99,...,1. nn110n由上式可知求In时,In1的误差的影响被缩小了。n=100时I100的近似值为0。

matlab 代码为

fprintf('稳定算法:\\n') y0=0; n=100;

plot(n,y0,'r*'); hold on

fprintf('y[100]=%10.6f',y0); while(1)

y1=1/10*[(1-exp(-n))/n-y0];

fprintf('y[%10.0f]=%10.6f',n-1,y1);plot(n-1,y1,'r*') if(n<=1) break;end y0=y1;n=n-1;

if mod(n,3)==0,fprintf('\\n'),end,end (具体值已省略) 编程实现得下图。 0.060.050.040.030.020.0100102030405060708090100

由图可知,该算法是稳定的。

1.3绘制静态和动态的Koch分形曲线

Koch曲线程序koch.m

function koch(a1,b1,a2,b2,n)

%koch(0,0,9,0,3)

%a1,b1,a2,b2为初始线段两端点坐标,n为迭代次数 %例如a1=0;b1=0;a2=9;b2=0;n=3;

%第i-1次迭代时由各条线段产生的新四条线段的五点横、纵坐标存储在数组A、B中 [A,B]=sub_koch1(a1,b1,a2,b2); for i=1:n

for j=1:length(A)/5;

w=sub_koch2(A(1+5*(j-1):5*j),B(1+5*(j-1):5*j)); for k=1:4

[AA(5*4*(j-1)+5*(k-1)+1:5*4*(j-1)+5*(k-1)+5),BB(5*4*(j-1)+5*(k-1)+1:5*4*(j-1)+5*(k-1)+5)]=sub_koch1(w(k,1),w(k,2),w(k,3),w(k,4)); end end A=AA; B=BB; end

plot(A,B) hold on axis equal

%由以(ax,ay),(bx,by)为端点的线段生成新的中间三点坐标并把这五点横、纵坐标依次分别存%储在数组A,B中

function [A,B]=sub_koch1(ax,ay,bx,by) cx=ax+(bx-ax)/3; cy=ay+(by-ay)/3; ex=bx-(bx-ax)/3; ey=by-(by-ay)/3;

L=sqrt((ex-cx).^2+(ey-cy).^2); alpha=atan((ey-cy)./(ex-cx)); if (ex-cx)<0

alpha=alpha+pi; end

dx=cx+cos(alpha+pi/3)*L; dy=cy+sin(alpha+pi/3)*L; A=[ax,cx,dx,ex,bx]; B=[ay,cy,dy,ey,by];

%把由函数sub_koch1生成的五点横、纵坐标A,B顺次划分为四组,分别对应四条折线段中 %每条线段两端点的坐标,并依次分别存储在4*4阶矩阵k中,k中第i(i=1,2,3,4)行数字代表第%i条线段两端点的坐标 function w=sub_koch2(A,B) a11=A(1);b11=B(1); a12=A(2);b12=B(2);

a21=A(2);b21=B(2); a22=A(3);b22=B(3); a31=A(3);b31=B(3); a32=A(4);b32=B(4); a41=A(4);b41=B(4); a42=A(5);b42=B(5);

w=[a11,b11,a12,b12;a21,b21,a22,b22;a31,b31,a32,b32;a41,b41,a42,b42];

%调用函数得到下图

n=5;

i=0;while ifigure(i+1);koch(0,0,3,0,i) i=i+1;pause(1)

end

1.5 1.5110.50.500-0.5-0.500.511.522.5300.511.522.531.51.5110.50.500-0.5-0.500.511.521.52.5300.511.522.5310.50

-0.500.511.522.53将pause(1)去掉可得静态图

2.1小行星轨道问题

为了确定方程中的5个待定系数,需要将上述5个点的坐标代入上面的方程

a1x22a2xya3y22a4x2a5y1,得:

a1x122a2x1y1a3y122a4x12a5y1 122a1x22a2x2y2a3y22a4x22a5y2 122ax2axyay233332a4x32a5y3 1 13ax22axyay22ax2ay 124434445414a1x522a2x5y5a3y522a4x52a5y5 1将这一包含5个未知数的线性方程组,写成矩阵的形式

x122x2x322x4x252x1y12x2y22x3y32x4y42x5y5y12y22y32y42y522x12x22x32x42x52y1a12y2a22y3a3=

2y4a42y5a511111AXb

(1) 求解这一线性方程组,即可得到曲线方程的系数

Y0=[6062 11179 16954 23492 68894];

A=zeros(5);X0(1); for i=1:5

A(i,1)=X0(i)^2;A(i,2)=2*X0(i)*Y0(i);A(i,3)=Y0(i)^2;A(i,4)=2*X0(i

);A(i,5)=2*Y0(i); end;

format long g;A

A=

2873496025 649907020 36747844 107210 12124 3417571600 1307048680 124970041 116920 2235 3951253881 2131422972 287438116 125718 33908 4443822244 3132047408 551874064 133324 46984 4746383236 9492766472 4746383236 137788 137788

B=[-1 -1 -1 -1 -1]';format long g;x=A\\B

x =

X0=[53605 58460 62859 66662 68894];

-8.06820280371841e-011 -7.63620099622306e-011 -3.0801152978055e-010 -8.89025159419867e-006 2.02829368401655e-005

(2)用Lu分解法解可得

format long g;

A = [2873496025 649907020 36747844 107210 12124; 3417571600 1307048680 124970041 116920 22358; 3951253881 2131422972 287438116 125718 33908; 4443822244 3132047408 551874064 133324 46984; 4746383236 9492766472 4746383236 137788 137788];

B=[-1 -1 -1 -1 -1]';

[L,U,flag ]=LU_Decom(A),format long g;x=U\\(L\\B)

flag = OK x =

-8.06820280370254e-011 -7.63620099622621e-011 -3.08011529780421e-010 -8.89025159420231e-006 2.02829368401615e-005 jacobi迭代法: jacobic(A)

因为谱半径不小于1,所以Jacobi迭代序列发散,谱半径SRH和B的所有特征值H如下: SRH =

4.41963931714337 ans =

-4.41963931714337 1.5453292801696 0.757138763648732 1.11838895650533 0.9987823168197

GSC(A)

因为谱半径不小于1,所以Gauss-Seidel迭代序列发散,谱半径SRH和B的所有特征值H如下: SRH =

1.12218280703645 ans =

0 0.0914047045360394 1.10703104191813 + 0.183783907450762i 1.10703104191813 - 0.183783907450762i 0.99748223605848

2.2

(1) 用Gauss列主元消去法、Gauss按比例列主元消去法、Cholesky分解求解下列线性方程组,并彼此互相验证。

(2) 判断用Jacobi 迭代法、Gauss-Seidel迭代法、SOR法(分别取. 若收敛,再用Jacobi 迭代法、0.8,1.2,1).解下列线性方程组的收敛性3Gauss-Seidel迭代法、SOR法(分别取0.8,1.2,1.3,1.6)分别解线性方程组,并比较各种方法的收敛速度.

x1x22x3x4x13x23x42x19x36x4x13x26x319x4(3) 用Cholesky分解求解下列线性方程组 4222410221430200

1,3,5,7.

4021812243344121411681344115310163340x10x620256x32033x423

103x5914x622142x715x2194580(1)

A=[1 -1 2 1;-1 3 0 -3;2 0 9 -6;1 -3 -6 19];b=[1 3 5 7]';

[RA,RB,n,x]=liezy(A,b), [RA,RB,n,x]=bilizy(A,b), cholesky(A,b) 列主元

因为RA=RB=n,所以此方程组有唯一解. RA = 4 RB = 4 n = 4 x =

-8.000000000000005 0.333333333333332 3.666666666666668 2.000000000000000 比例主元

因为RA=RB=n,所以此方程组有唯一解. RA = 4 RB = 4 n = 4 x =

-8.000000000000000 0.333333333333333 3.666666666666667 2.000000000000000 cholesky分解 S1 = 0 S1 = 0 S1 = 0 x =

0 0 3.666666666666673 2.000000000000003 x =

0 0.333333333333329 3.666666666666673 2.000000000000003 x =

-8.000000000000020 0.333333333333329 3.666666666666673 2.000000000000003

-8.000000000000020 0.333333333333329 3.666666666666673 2.000000000000003 对比可知,三种方法可相互验证。

(2)

x0=ones(4,1);[]

A=[1 -1 2 1;-1 3 0 -3;2 0 9 -6;1 -3 -6 19]; jacobic(A),GSC(A) 因为谱半径小于1,所以Jcobi迭代序列收敛,谱半径SRH和B的所有特征值H如下: SRH =

0.966881024532242 ans =

0.966881024532242 0.536020655954138 -0.903317605697383 -0.599584074788997

因为谱半径小于1,所以Gauss-Seidel迭代序列收敛,谱半径SRH和B的所有特征值H如下: SRH =

0.934888367800313 ans =

0 0.934888367800313 0.281485901205534

-2.38957160758552e-017 jacobi迭代法:

A=[1 -1 2 1;-1 3 0 -3;2 0 9 -6;1 -3 -6 19]; b=[1 3 5 7]';x0=ones(4,1); [x,n]=jacobi(A,b,x0) x =

-7.99997357113847 0.333339142428815 3.66665839089021 1.99999680708368 n = 376

gauseidel法

[x,n]=gauseidel (A,b,x0) x =

-7.99998673303378 0.333336112109258 3.66666262275452 1.99999846346783 n = 197

SOR法:

w1=0.8;w2=1.2;w3=1.3;w4=1.6;

SORC(A,w1), SORC(A,w2), SORC(A,w3), SORC(A,w4)

w =

0.8

因为谱半径小于1,所以SOR迭代序列收敛,谱半径SRH和B的所有特征值H如下: SRH =

0.956498743812875 ans =

0.956498743812875 0.503328665301619 0.050190208900126 0.0662163001140356 w =

1.2

因为谱半径小于1,所以SOR迭代序列收敛,谱半径SRH和B的所有特征值H如下: SRH =

0.901948033909749 ans =

0.901948033909749 0.0392206356524567 0.00773145469258105 + 0.212532205268192i 0.00773145469258105 - 0.212532205268192i w =

1.3

因为谱半径小于1,所以SOR迭代序列收敛,谱半径SRH和B的所有特征值H如下: SRH =

0.877539492310148 ans =

0.877539492310148 0.0946172161251248 -0.0537947284866427 + 0.307669991725664i -0.0537947284866427 - 0.307669991725664i w =

1.6

因为谱半径小于1,所以SOR迭代序列收敛,谱半径SRH和B的所有特征值H如下: SRH =

0.608878775828679 ans =

0.590094939083101 + 0.0369506212620504i 0.590094939083101 - 0.0369506212620504i -0.21966219054509 + 0.56787488560383i -0.21966219054509 - 0.56787488560383i w=0.8时 [x,n]=SOR(A,b,x0,0.8,10^(-5),300) x =

-7.99980126177984 0.333375649565234 3.66660553355675 1.99997665093436

n = 238

w=1.2时

[x,n]=SOR(A,b,x0,1.2,10^(-5),300) x =

-7.99992226324999 0.333349197733832 3.66664330811682 1.99999119661784 n = 111

w=1.3时

[x,n]=SOR(A,b,x0,1.3,10^(-5))

x =

-7.99993630689315 0.333346074182271 3.66664773570359 1.99999290958378 n = 89

w=1.6时

[x,n]=SOR(A,b,x0,1.6,10^(-5))

x =

-7.99999410032529 0.333332406655426 3.66666433369261 1.99999911550516 n = 26

比较上述方法可得,w=1.6时的SOR法最为快速。

(3)Cholesky分解求解下列线性方程组

A=[4 2 -4 0 2 4 0 0;2 2 -1 -2 1 3 2 0;-4 -1 14 1 -8 -3 5 6; 0 -2 1 6 -1 -4 -3 3;2 1 -8 -1 22 4 -10 -3;4 3 -3 -4 4 11 1 -4;

0 2 5 -3 -10 1 14 2;0 0 6 3 -3 -4 2 19];b=[0;-6;20;23;9;-22;-15;45];

Choleshy(A,b)

121.1481 -140.1127 29.7515 -60.1528 10.9120 -26.7963 5.4259 -2.0185

2.3设

2023100,b0,b1,b0 A1811232315001(1) 将矩阵A进行LU分解A=LU,其中U是上三角矩阵,L是主对角线上的元素都是1的下三角矩阵。

(2) 利用上述分解分别求解方程组

AXb1,1AXb2,AXb3

并由此求出逆矩阵A。

(3) 用LU分解求下列线性方程组的解

4238654220214268680211610114620010x15x12021x334x421673323x53

57172635x6463425301x713917342122x838713920124x9198324863121x10123613511523001101010039(方程组的精确解是x(1,1,0,1,2,0,3,1,1,2)T.)

(1)

A =[20 2 3;1 8 1;2 -3 15] ;[L,U,flag ]=LU_Decom(A)

b1=[1 0 0]'; b2=[0 1 0]'; b3=[0 0 1]';x1=A\\b1, x2=A\\b2, x3=A\\b3

L =

1.0000 0 0 0.0500 1.0000 0 0.1000 -0.4051 1.0000 U =

20.0000 2.0000 3.0000 0 7.9000 0.8500 0 0 15.0443 flag = OK

X=[x1,x2,x3] X =

0.0517 -0.0164 -0.0093 -0.0055 0.1237 -0.0072 -0.0080 0.0269 0.0665 X即为A的逆。

(2) matlab中输入命令:[L,U,flag ]=LU_Decom(A), format long g;x=U\\(L\\b)得

L为

U为

flag = OK x =

1 -0.999999999999999 4.21884749357559e-015 0.999999999999999

1.99999999999999 1.68753899743024e-015 3 0.999999999999999 -1 2

2.4 (1) 用追赶法求解方程组AXb

41(a) A14114121,b 11430302301521125210125210(b) A,b 125210125202011252020(2) 设计实验验证Hilbert矩阵Hn的病态性,其中

11Hn21n

12131n11n1n1. 12n1(1)

a1=4*ones(30,1)';a2=ones(29,1)';a3=a2;b=ones(30,1)';b(1,1)=2;b(30,1)=2; A=diag(a1)+diag(a3,1)+diag(a2,-1);

x =chasing (A, b)

zhuiganfa(A,b)

ans = 0.5359 -0.1436 0.0385

-0.0103 0.0028 -0.0007 0.0002 -0.0001 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0001 0.0002 -0.0007 0.0028 -0.0103 0.0385 -0.1436 0.5359

(b)

a1=5*ones(20,1)';a2=-2*ones(19,1)';a3=a2;a4=ones(18,1)';a5=a4;b=zeros(20,1);b(1,1)=1;

A=diag(a1)+diag(a3,1)+diag(a2,-1)+diag(a4,-2)+diag(a5,2);zhuiganfa(A,b) ans =

0.249999999999829 0.124999999999574 0.0624999999991047 0.0312499999981881 0.0156249999963656 0.00781249999272582 0.00390624998544897 0.00195312497089661 0.000976562441792561 0.000488281133584789

0.000244140392169412 0.00012206984683874 6.10342249274393e-005 3.05157154798577e-005 1.5255063772205e-005 7.62194395065481e-006 3.79979610443202e-006 1.87754631042523e-006 8.94069671631063e-007 3.57627868652425e-007

(2)

matlab代码为:

m=input('input m:='); %输入矩阵的阶数 N=[m]; N=[m];

for k=1:length(N) n=N(k); %矩阵的阶 H=hilb(n); %产生n阶Hilbert矩阵

disp(H) %输出n阶Hilbert矩阵 Hi=invhilb(n); %产生完全准确的n阶逆Hilbert矩阵 b=ones(n,1); %生成n阶全1向量

x_approx=H\\b; %利用左除H求近似解

x_exact=Hi*b; %利用准确逆Hilbert矩阵求准确解

ndb=norm(H*x_approx-b);nb=norm(b); ndx=norm(x_approx - x_exact);

nx=norm(x_approx);

er_actual(k)=ndx/nx; %实际相对误差

K=cond(H); %计算Hilbert矩阵的条件数 er_approx(k)=K*eps; %最大可能的近似相对误差 er_max(k)=K*ndb/nb; %最大可能的相对误差 end

disp('Hilbert矩阵阶数'),disp(N), format short e disp('实际误差 er_actual'),disp(er_actual),disp('')

disp('近似的最大可能差 er_approx'),disp(er_approx),disp('') disp('最大可能误差 er_max'),disp(er_max),disp('')

执行代码,在matlab中输入4得 input m:=

1.0000e+000 5.0000e-001 3.3333e-001 2.5000e-001 5.0000e-001 3.3333e-001 2.5000e-001 2.0000e-001 3.3333e-001 2.5000e-001 2.0000e-001 1.6667e-001 2.5000e-001 2.0000e-001 1.6667e-001 1.4286e-001 Hilbert矩阵阶数 4

实际误差 er_actual 1.0284e-013

近似的最大可能差 er_approx 3.4447e-012 最大可能误差 er_max 4.7732e-011

多次执行代码,分别将m值设为8 10 16 100后得

阶数 实际误差 近似的4 8 10 16 100 1.0284e-013 1.7310e-007 1.9489e-004 2.6458e+002 9.1294e+123 最大可3.4447e-012 3.3879e-006 3.5583e-003 能误差 最大可能误差 1.3948e+002 3.5120e+004 4.7732e-011 3.8709e-002 1.2703e+003 6.2587e+009 1.2236e+013 对于高阶系数矩阵来说,如果系数矩阵是Hilbert矩阵,则迭代结果误差较大,而且阶数越大,误差就越大。

3.1

用Taylor级数的第i项多项式pi(x)(i1,2,图形。

,30),分别在区间[0,2]和

[0,]上逼近正弦函数f(x)sinx,x[0,],并用计算机绘出上面31个函数的

syms x; y=sin(x);

ezplot(y,[0,pi/2]); grid on;

axis([0,pi/2,0,3]); hold on; for m=1:2:30

p=taylor(y,x,m); ezplot(p,[0,pi/2]); axis([0,pi/2,0,1.5]); end

xlabel('x'); ylabel('y'); grid on

- x27/10888869450418352160768000000 +...+ x1.51y0.5000.5x11.5

3.2追赶曲线的计算机模拟

算法分析

(1)当狼位于点Pk(xk,yk)时,兔子位于点Qk(0,0tk),则狼在[tk,tk1]时段内的运动方向为PkQk(xk,v0tkyk),所以ekPkQk/norm(PkQk)

(2)PkPk1eknorm(PkPk1)2ekv0(tk1tk),(xk1xk,yk1yk)2v0(tk1tk)ek

程序及运行结果

Q=[0 0];%兔子坐标 P=[100 0];%狼坐标 PQ=Q-P;%狼的方向向量

step =1;%模拟步长:兔子奔跑的距离,step越小就越精确 count = 60/step;%以兔子的奔跑距离划分 PQ=PQ/norm(PQ)*step;%归一化,单位向量 trackP=P; trackQ=Q;

for k=1:count; P = P + 2*PQ;%2倍速度

Q = Q + step*[0 1];%[0 1]为兔子奔跑方向的单位方向向量 PQ = Q - P;

trackP(1+k,:)=P; trackQ(1+k,:)=Q;

PQ=PQ/norm(PQ)*step;%归一化,单位向量 dis=sqrt(sum((P-Q).^2));

plot(trackP(:,1),trackP(:,2),Q(1),Q(2),'rp',0,60,'r+'); pause(0.02) end

dis%兔子到达窝边时,狼兔之间的距离 P %兔子到达窝边时,狼的坐标 Q %兔子到达窝边时,兔子的坐标

运行结果: dis =

7.061909001184961 P =

1.680502483487450 53.140957053348501 Q =0 60

(将pause(0.2)去掉即为静态模型)

60504030201000102030405060708090100

4.1方程xsinx1的根实际上是两个函数y1sinx,y21x的交点的横坐标,用

计算机绘出两个函数在区间[6,6]的图形,观察图形,分析它们的交点分布规律

及特点,试写出方程的全部实根所在的隔根区间;并用二分法求出每一个根的近似值。

定义二分法代码:eff.m

function eff(fun,a,b,eps) if nargin<4 eps=1e-5; end

fa=feval(fun,a);fb=feval(fun,b); if fa*fb>0

disp('[a,b]不是有根区间,请重新调整'); return end

while abs(b-a)/2>eps

x=(a+b)/2;fx=feval(fun,x); if fx*fa<0 b=x;fb=fx; else

a=x;fa=fx; end end

x=(a+b)/2 主程序:

x=-6:0.1:6;

y1=sin(x);y2=1./x;plot(x,y1,'b',x,y2,'r');

151050-5-10-15-6-4-20246

通过函数图像可知,在0,上y1/x为单调递减函数,而ysinx为周期函数 第一个交点在0,/2之间,之后的方程的全部实根所在的隔根区间为

k/2,k3/2

调用eff.m

fun=inline('x*sin(x)-1');eff(fun,0,pi/2) for k=0:1

eff(fun,k*pi+pi/2,k*pi+3*pi/2) end ans = 1.1141 ans = 2.7725 ans =

6.4391

4.2设方程f(x)x33x10有3个实根

*x20.34727,*x31.53209

*x11.8793,**尽可能多地采用下面六种不同计算格式,求f(x)0的根x1*或x2或x3的近似值,

并观察相应格式的收敛性。

3x1(1) xk1k2,k0,1,2,xk3xk1; (2) xk1,k0,1,2,3;

(3) xk13xk1,k0,1,2,3; (4) xk11,k0,1,2,2xk3;

1(5) xk13,k0,1,2,xk33xk11xk; (6) xk1xk,k0,1,2,. 23xk1定义函数ddf

function [x,k]=ddf(fun,x0,eps,Nmax) if nargin<4 Nmax=500; end

if nargin<3 eps=1e-5; end

x=x0;x0=x+2*eps;k=0; while abs(x0-x)>eps&kif k==Nmax

warning('已迭代次数上限 '); end x

运行结果:

(1)fun=inline('(3*x+1)/x^2');[x,k]=ddf(fun,1.8,1e-5);

x = NaN

在1.8发散

fun=inline('(3*x+1)/x^2');[x,k]=ddf(fun,-0.3,1e-5)

x = NaN x = NaN k =

31

在-0.3发散

fun=inline('(3*x+1)/x^2');[x,k]=ddf(fun,-1.5,1e-5) x =

-1.5321 x =

-1.5321 k =

28

在-1.5收敛 (2)

fun=inline('(x^3-1)/3');[x,k]=ddf(fun,1.8,1e-5) x = -0.3473 x =-0.3473 k =9

在1.8收敛

fun=inline('(x^3-1)/3');[x,k]=ddf(fun,-0.3,1e-5) x =-0.3473 x =-0.3473 k = 5

在-0.3收敛

fun=inline('(x^3-1)/3');[x,k]=ddf(fun,-1.5,1e-5) x =-0.3473 x =-0.3473 k = 12

在-1.5收敛

(3) fun=inline('(3*x+1)^(1/3)');[x,k]=ddf(fun,1.8,1e-5)

x =1.8794 x =1.8794 k = 8

在1.8收敛

fun=inline('(3*x+1)^(1/3)');[x,k]=ddf(fun,-0.3,1e-5) x = 1.8794 x =1.8794 k =12

在-0.3收敛

fun=inline('(3*x+1)^(1/3)');[x,k]=ddf(fun,-1.5,1e-5) x =1.8794 + 0.0000i x =1.8794 + 0.0000i k =12

在-1.5收敛

(4)fun=inline('1/(x^2-3)');[x,k]=ddf(fun,1.8,1e-5)

x =-0.3473 x =-0.3473 k =7

在1.8收敛

fun=inline('1/(x^2-3)');[x,k]=ddf(fun,-0.3,1e-5) x =-0.3473 x =-0.3473 k = 5

在-0.3收敛

fun=inline('1/(x^2-3)');[x,k]=ddf(fun,-1.5,1e-5)

x =-0.3473 x =-0.3473 k = 8

在-1.5收敛 (5)

fun=inline('(3+1/x)^(1/3)');[x,k]=ddf(fun,1.8,1e-5)

x = 1.5396 x =1.5396 k = 5

在1.8收敛

fun=inline('(3+1/x)^(1/3)');[x,k]=ddf(fun,-0.3,1e-5) x = 1.5396 + 0.0000i x =1.5396 + 0.0000i k = 7

在-0.3收敛

fun=inline('(3+1/x)^(1/3)');[x,k]=ddf(fun,-1.5,1e-5) x = 1.5396 x =1.5396 k = 6

在-1.5收敛 (6)

fun=inline('x-(1/3)*(x^3-3*x-1)/(x^2-1)');[x,k]=ddf(fun,-1.5,1e-5) x =-1.5321 x =-1.5321 k = 3

在-1.5收敛

fun=inline('x-(1/3)*(x^3-3*x-1)/(x^2-1)');[x,k]=ddf(fun,-0.3,1e-5) x = -0.3473 x =-0.3473 k = 3

在-0.3收敛

fun=inline('x-(1/3)*(x^3-3*x-1)/(x^2-1)');[x,k]=ddf(fun,1.8,1e-5) x = 1.8794 x =1.8794 k= 4

在1.8收敛

4.3一个10次多项式 P(x)x[1,a1,a2,10a1x9a2x8a9xa10的系数为

,a9,a10]=[1, -55, 1320, –18150, 157773, –902055, 3416930,

–8409500, 12753576, –10628640, 6328800]

(1) 用多项式的求根指令roots求出方程P(x)0的10个根;

(2) 试设计一个算法,求出方程P(x)0的10个根的近似值,并与(1)的结果进行比较;

(3) 将方程P(x)0左边多项式P(x)中的9次项的系数–55改为–56,得到一个新的10次方程;求解新的方程,观察根的变化是否很显著,说明了什么?

(1)

roots([1 -55 1320 -18150 157773 -902055 3416930 -8409500 12753576 -10628640 6328800]) ans =

10.6051 + 1.0127i 10.6051 - 1.0127i 8.5850 + 2.7898i 8.5850 - 2.7898i 5.5000 + 3.5058i 5.5000 - 3.5058i 2.4150 + 2.7898i 2.4150 - 2.7898i 0.3949 + 1.0127i 0.3949 - 1.0127i

(2)

原代码为:

function root=Parabola(f,a,b,x,eps) if(nargin==4) eps=1.0e-4; end

f1=subs(sym(f),findsym(sym(f)),a); f2=subs(sym(f),findsym(sym(f)),b); if(f1==0) root=a; end if(f2==0) root=b; end

if(f1*f2>0)

disp('Á½¶Ëµãº¯ÊýÖµ³Ë»ý´óÓÚ0!'); return; else tol=1;

fa=subs(sym(f),findsym(sym(f)),a); fb=subs(sym(f),findsym(sym(f)),b); fx=subs(sym(f),findsym(sym(f)),x); d1=(fb-fa)/(b-a); d2=(fx-fb)/(x-b); d3=(d2-d1)/(x-a); B=d2+d3*(x-b);

root=x-2*fx/(B+sign(B)*sqrt(B^2-4*fx*d3)); t=zeros(3); t(1)=a; t(2)=b; t(3)=x; while(tol>eps) t(1)=t(2); t(2)=t(3); t(3)=root;

f1=subs(sym(f),findsym(sym(f)),t(1)); f2=subs(sym(f),findsym(sym(f)),t(2)); f3=subs(sym(f),findsym(sym(f)),t(3)); d1=(f2-f1)/(t(2)-t(1)); d2=(f3-f2)/(t(3)-t(2)); d3=(d2-d1)/(t(3)-t(1)); B=d2+d3*(t(3)-t(2));

root=t(3)-2*f3/(B+sign(B)*sqrt(B^2-4*f3*d3)); tol=abs(root-t(3)); end end

f=inline(’

x^10-55*x^9+1320*x^8-18150*x^7+157773*x^6-902055*x^5+3416930*x^4-8409500*x^3+12753576*x^2-10628640*x^1+6328800‘);

Parabola('

x^10-55*x^9+1320*x^8-18150*x^7+157773*x^6-902055*x^5+3416930*x^4- 8409500*x^3+12753576*x^2-10628640*x+6328800', 10.6051, 10.6051 - 1.0128i,1.0e-8) ans =

2.41497968303632 + 2.7898128565398i

(3)

roots([1 -56 1320 -18150 157773 -902055 3416930 -8409500 12753576 -10628640 6328800]) ans =

21.7335121408968 7.35013402487903 + 7.79728509446717i 7.35013402487903 - 7.79728509446717i 4.35886185227404 + 3.3285467123814i 4.35886185227404 - 3.3285467123814i 5.1830729717618 2.4378000685152 + 2.79739872415432i 2.4378000685152 - 2.79739872415432i 0.394911498002447 + 1.01269497774644i

0.394911498002447 - 1.01269497774644i 原多项式方程的根 10.6051 + 1.0127i 10.6051 - 1.0127i 8.5850 + 2.7898i 8.5850 - 2.7898i 5.5000 + 3.5058i 5.5000 - 3.5058i 2.4150 + 2.7898i 2.4150 - 2.7898i 0.3949 + 1.0127i 0.3949 - 1.0127i 变化之后多项式方程的根 21.7335 7.3501 + 7.7973i 7.3501 - 7.7973i 4.3589 + 3.3285i 4.3589 - 3.3285i 5.1831 2.4378 + 2.7974i 2.4378 - 2.7974i 0.3949 + 1.0127i 0.3949 - 1.0127i 将方程P(x)0左边多项式P(x)中的9次项的系数–55改为–56,得到一个新的10次方程;求解新的方程得根的变化较为显著。说明次数越高的项的 微小变化能够引起根的巨大变化。

4.4找出方程x32x50的有根区间,再分别用二分法、不动点迭代法(取迭

代函数1(x)352x)、Steffensen迭代法、Newton迭代法、弦截法分别求方程

x32x50的根, 要求|xk1xk|106.

fplot('[x^3+2*x-5,0]',[1,2]);grid

76543210-1-211.11.21.31.41.51.61.71.81.92

由上图可知方程根的区间为(1.3,1.4).

二分法:

root=HalfInterval('x^3+2*x-5',1.3,1.4) root = 1.3282470703125 不动点迭代:

先用M-文件写一个名为g.m的函数 function y=g(x); y=(5-2*x)^(1/3); 调用fixpth函数

fixpt('g',1.3,10^(-6),30)

k =2 err =0.0389 k =3 err =0.0146 k =4 err =0.0055 k =5 err =0.0021 k =6 err =7.8957e-004 k =7 err =2.9832e-004 k =8 err = 1.1273e-004 k =9 err =4.2596e-005 k =10 err =1.6095e-005 k =11 err =6.0819e-006 k =12 err =2.2981e-006 k =13 err =8.6839e-007 P =

1.3000 1.3389 1.3243 1.3298 1.3277 1.3285 1.3282 1.3283 1.3283 1.3283 1.3283 1.3283 1.3283 ans =

1.3000

Steffensen迭代法:

[x_steff,time_steff]=Steff(' x^3+2*x-5',1.3) x_steff = 1.3283 time_steff = 4 newton法:

用一个名为f.m的文件定义函数 f(x)=x^3+2x-5; 后用一个名为 df.m的文件定义函数的微商 df=3x^2+2

newton('f','df',1.3,10^(-6),20) p0 =1.3 ans = -0.202999999999999

p1 = 1.32871287128713 err = 0.0287128712871285 k =1 y = 0.00323894473556585 p1 = 1.3282689633622 err =0.000443907924925657 k =2 y =7.85398064806486e-007 p1 =1.32826885566861 err =1.07693588269342e-007 k =3

y =4.52970994047064e-014 ans =1.32826885566861

弦截法:

secant('f',1.2,1.3,10^(-6),20 )

k = 0 p0 = 1.2 p1 = 1.3 ans =-0.872

ans = -0.202999999999999 k =1

p1 =1.33034379671151 err =0.0303437967115097 y =0.0151494910753769 k =2

p1 = 1.32823655788547 err =0.00210723882603614 y =-0.000235540166951154 k = 3

p1 =1.32826881907345

err =3.22611879737256e-005 y =-2.66884642385889e-007 k =4

p1 =1.32826885566925 err =3.65958068293537e-008 y =4.70912198125006e-012 ans =1.32826885566925

5.1

x

0.4 0.55 0.65 0.8 0.9 1.05 0.41075 0.57815 0.69675 0.88811 1.02652 1.25382

f(x)

(1) 分别用Lagrange插值和Newton插值求f(0.596)的近似值; (2) 分别用Aitken算法和Neville算法求f(0.596)的近似值;

(3) 分别用Newton向前插值公式和Newton向后插值公式求f(0.42)和

f(1.0)的近似值。

(1)

lagrange法

format long

g,x=[0.4,0.55,0.65,0.8,0.9,1.05];y=[0.41075,0.57815,0.69675,0.88811,1.02652,1.25382];

Language(x, y),Language(x, y,0.596)

ans =

0.00029304*t^5 + 0.0302711*t^4 + 0.123615*t^3 + 0.0296166*t^2 + 0.990118*t + 0.0012748 ans =

0.631917499231745

newton插值法

f=Newtonchazhi(x,y),a=Newtonchazhi(x,y,0.596)

f =

0.00029304*t^5 + 0.0302711*t^4 + 0.123615*t^3 + 0.0296166*t^2 + 0.990118*t + 0.0012748 a =

0.631917499231745

(2)

Atken 法

f=Atken(x,y),a=Atken(x,y,0.596)

f =

0.00029304*t^5 + 0.0302711*t^4 + 0.123615*t^3 + 0.0296166*t^2 + 0.990118*t + 0.0012748 a =

0.631917499231746

neville法:

neville(x,y,0.596)

ans =

0.631917499231746

(3)

newton向前差分

f=Newtonforward(x,y),a= Newtonforward(x,y,0.42),b=

Newtonforward(x,y,1.0) f =

0.00429017*t^5 - 0.0532046*t^4 + 0.232233*t^3 - 0.41302*t^2 + 0.397101*t + 0.41075 a =

0.520282035525315 b =

0.57815

newton向后差分

a= Newtonbackward(x,y,0.42), b=Newtonbackward(x,y,1.0)

a =

-25.36037050752 b =

1.14756289437586

5.2

已知函数f(x)cosx在节点x0,0.25,0.5,0.75,1.0处的值,求 (1) 满足边界条件f(0)0和f(1.0)0的三次样条插值函数S(x); (2) 计算S(x)dx,并计算f(x)dxS(x)dx;

000111(3) 用S(0.5)和S(0.5)分别近似f(0.5)和f(0.5).

(1)

x=[0,0.25,0.5,0.75,1.0] ;y=cos(pi.*x);dx0=0;dxn=0;s=csfit(x,y,dx0,dxn) s =

2.0281 -5.1933 0 1.0000 4.8963 -3.6722 -2.2164 0.7071 4.8963 -0.0000 -3.1344 0.0000

2.0281 3.6722 -2.2164 -0.7071

即求得:

s12.0281x35.1933x21 x[0,0.25]

s24.8963x33.6722x22.2164x0.7071x[0.25,0.5] s34.8963x33.1344x[0.5,0.75]

s42.0281x33.6722x22.2164x0.7071x[0.5,0.75]

(2)

syms x;

y1=2.0281*x^3-5.1933*x^2+1; int(y1,0,0.25) ans =

460661/2048000 a1=vpa(ans) a1 =

-0.178949707031255

y2=4.8963*x^3-3.6722*x^2-2.2164*x+0.7071; int(y2,0.25,0.5) ans =

-2862233/30720000 a2=vpa(ans) a2 =

-0.093171647135416666666666666666667

y3=4.8963*x^3-3.1344*x;

int(y3,0.5,0.75) ans =

-366489/2048000 a3=vpa(ans) a3 =

-0.17894970703125

y4=2.0281*x^3+3.6722*x^2-2.2164*x-0.7071; int(y4,0.75,1) ans =

12062213/30720000 a4=vpa(ans)

a4 =

0.39265016276041666666666666666667 a1+a2+a3+a4 ans =

-0.0584208984375 即为

10S(x)dx

10f(x)dx为

int(cos(pi*x),0,1) ans = 0

101f(x)dxS(x)dx=0.0584208984375

0(3)

diff(4.8963*x^3-3.6722*x^2-2.2164*x+0.7071) ans =

(146889*x^2)/10000 - (18361*x)/2500 - 5541/2500

5.3

xi yi yi 0 1 2 3 4 5 6 7 8 9 10 0.0 0.79 1.53 2.19 2.71 3.03 3.27 2.89 3.06 3.19 3.29 0.8 0.2 (1)试求出三次样条插值曲线, 并画出图形。

(2) 数据中的一阶导数的信息换成s(0)s(10)0,试求出三次样条插值曲线, 并画出图形。

(1)

调用csfit函数

x=0:10;

y=[0.0 0.79 1.53 2.19 2.71 3.03 3.27 2.89 3.06 3.19 3.29]; dx0=0.8 ;dxn=0.2;csfit(x,y,dx0,dxn) ans =

-0.00851403685003569 -0.00148596314996432 0.8 0 -0.00445788944989297 -0.0270280737000714 0.771485963149964 0.79 -0.00365440535039247 -0.0404017420497503 0.704056147400143 1.53 -0.040924489148537 -0.0513649581009277 0.612289447249465 2.19 0.10735236194454 -0.174138425546539 0.386786063601998 2.71 -0.268484958629623 0.147918660287082 0.360566298342541 3.03 0.426587472573951 -0.657536215601787 -0.149051256972164 3.27 -0.267864931666182 0.622226202120067 -0.184361270453885 2.89 0.054872254090777 -0.18136859287848 0.256496338787702 3.06 0.0583759153030742 -0.0167518306061485 0.0583759153030744 3.19

函数图像为:

543210-1-2012345678910

(2)

x=0:10;

y=[0.0 0.79 1.53 2.19 2.71 3.03 3.27 2.89 3.06 3.19 3.29]; yangtiao2(x,y) 第二种边界条件:

543210-1-2012345678910

5.4给定初值问题

y2yt2et,1t2,ty(1)0,

其精确解为y(t)t2ete:

(1) 取h0.1,用4阶Runge-Kutta方法求上述初值问题的近似解,并将它们与y的实际值进行比较;

(2) 用(1)中所得的答案分别用分段线性插值、分段三次Hermite插值、三次样条求在下列值处的近似y,并将它们与y的实际值进行比较。

( i ) y1.04, (ii) y1.55, (ⅲ) y1.97.

(1)

四阶龙格库塔公式 x=[];y=[];

a=1;b=2;h=0.1;alpha=0;f=inline('y*2/t+t^2*exp(t)'); x=a:h:b; n=(b-a)/h; y(1)=alpha;

disp('四阶龙格-库塔法求得的解为:') for i=2:n+1

k1=feval(f,x(i-1),y(i-1));

k2=feval(f,x(i-1)+h/2,y(i-1)+k1*h/2); k3=feval(f,x(i-1)+h/2,y(i-1)+k2*h/2); k4=feval(f,x(i-1)+h,y(i-1)+k3*h);

y(i)=y(i-1)+h/6*(k1+2*k2+2*k3+k4); plot(x(i),y(i),'+') ,hold on;

sprintf('x=%3.1f,y=%9.7f',x(i),y(i)) end

四阶龙格-库塔法求得的解为: ans =

x=1.1,y=0.3459103 ans =

x=1.2,y=0.8666217 ans =

x=1.3,y=1.6071813 ans =

x=1.4,y=2.6203113 ans =

x=1.5,y=3.9676019 ans =

x=1.6,y=5.7208793 ans =

x=1.7,y=7.9637718 ans =

x=1.8,y=10.7935018 ans =

x=1.9,y=14.3229357 ans =

x=2.0,y=18.6829266

实际的值分别为:

x=1:0.1:2;y=x.^2.*(exp(x)-exp(1));plot(x,y)

2018161412108642011.11.21.31.41.51.61.71.81.92

‘+’表示用龙格-库塔法所求的解,曲线为真实解,由图可知龙格库塔法非常精确。

(2)

分段线性插值

x=1:0.1:2;

y=[0 0.3459103 0.8666217 1.6071813 2.6203113 3.9676019 5.7208793

7.9637718 10.7935018 14.3229357 18.6829266];

xi1=[1.04 1.55

1.97];yi1=interp1(x,y,xi1,'linear'),plot(xi1,yi1),hold on; yi1 =

0.1384 4.8442 17.3749

分段三次样条插值

yi2=interp1(x,y,xi1,'spline') yi2 =

0.1202 4.7885 17.2797

分段三次Hermite插值

y1= y.*2./x+x.^2.*exp(x),yi3=Hermite(x,y,y1,xi1) y1 =

2.7183 4.2640 6.2253 8.6737 11.6915 15.3739 19.8309 25.1889 31.5936 39.2129 48.2392 yi3 =

0.1159 4.7886 17.2838

实际值:

y=xi1.^2.*(exp(xi1)-exp(1)) y =

0.1200 4.7886 17.2793

对比可知样条插值更为精确

6.2设f(x)sin(x),求次数分别为2、3、6、8、10的多项式p(x)使得

10f(x)p(x)dxmin并且画出p(x)和f(x)的图像进行比较。

2

x=0:0.01:1;y=sin(pi.*x);z1=polyfit(x,y,2);

z2=polyfit(x,y,3);z3=polyfit(x,y,6);z4=polyfit(x,y,8);z5=polyfit(x,y,10);pz1=poly2str(z1,'x'), pz2=poly2str(z2,'x'),

pz3=poly2str(z3,'x'), pz4=poly2str(z4,'x'), pz5=poly2str(z5,'x'),

pz1 =

-4.1076 x^2 + 4.1076 x - 0.047495 pz2 =

-1.1085e-014 x^3 - 4.1076 x^2 + 4.1076 x - 0.047495 pz3 =

-1.2275 x^6 + 3.6825 x^5 - 0.55998 x^4 - 5.0175 x^3 - 0.020227 x^2 + 3.1427 x - 1.2757e-005 pz4 =

0.22025 x^8 - 0.88098 x^7 + 0.20958 x^6 + 2.4547 x^5 + 0.026512 x^4

- 5.172 x^3 + 0.00035819 x^2 + 3.1416 x + 8.2738e-008 pz5 =

-0.024436 x^10 + 0.12218 x^9 - 0.03994 x^8 - 0.57331 x^7 - 0.011209 x^6

+ 2.5534 x^5 - 0.00057975 x^4 - 5.1676 x^3 - 3.5932e-006 x^2 + 3.1416 x

- 3.4178e-010 ;

令:

y1= -4.1076.* x.^2 + 4.1076.*x - 0.047495;

y2=-1.1085e-014.* x.^3 - 4.1076.*x.^2 + 4.1076.*x - 0.047495;

y3=-1.2275.*x.^6 + 3.6825.*x.^5 - 0.55998.*x.^4 - 5.0175.*x.^3 - 0.020227.*x.^2+ 3.1427.*x - 1.2757e-005;

y4 = 0.22025.*x.^8 - 0.88098.*x.^7 + 0.20958.*x.^6 + 2.4547.*x.^5 + 0.026512.*x.^4 - 5.172.*x.^3 + 0.00035819.*x.^2 + 3.1416.*x + 8.2738e-008;

y5 =-0.024436.*x.^10 + 0.12218.*x.^9 - 0.03994.*x.^8 - 0.57331.*x.^7 - 0.011209.*x.^6+ 2.5534.*x.^5 - 0.00057975.*x.^4 - 5.1676.*x.^3 - 3.5932e-006.*x.^2 + 3.1416.*x- 3.4178e-010;

plot(x,y,x,y1,x,y2,x,y3,x,y4,x,y5); legend('原函数曲线','二次多项式拟合','三次多项式拟合','六次多项式拟合','八次多项式拟合','十次多项式拟合')

1.2原函数曲线二次多项式拟合三次多项式拟合六次多项式拟合八次多项式拟合十次多项式拟合 10.80.60.40.20-0.2 00.10.20.30.40.50.60.70.80.91

由上图可知,二次,三次多项式拟合的效果有一定误差,而三次,六次,八次,十次多项式拟合曲线与原曲线十分吻合,几乎重合,表明次数越高,拟合效果越好。

6.3神经元模型用于蠓的分类识别

画出15对数据的散点图,其中,Af用*标记,Apf用+标记。

x1=[1.24,1.36,1.38,1.38,1.38,1.40,1.48,1.54,1.56];y1=[1.27,1.74,1.64,1.82,1.90,1.70,1.82,1.82,2.08];

x2=[1.14,1.18,1.20,1.26,1.28,1.30];y2=[1.78,1.96,1.86,2.00,2.00,1.96];

plot(x1,y1,'*',x2,y2,'+')

2.121.91.81.71.61.51.41.31.11.151.21.251.31.351.41.451.51.551.6

分别求出两个超定方程组的系数:

x1=[1.24,1.36,1.38,1.38,1.38,1.40,1.48,1.54,1.56];y1=[1.27,1.74,1.64,1.82,1.90,1.70,1.82,1.82,2.08];

x2=[1.14,1.18,1.20,1.26,1.28,1.30];y2=[1.78,1.96,1.86,2.00,2.00,1.96];zxec(x1,y1,1);zxec(x2,y2,1)

alpha =

-0.8174 1.8197

alpha =

0.5757

1.1014

11.81972根据这两组数据可解得:w00.1735,w21.4356。由于我们取其平

11.10142均值,则

11.4606,12.0968,所以求得的方程为:22.0968x1.4356y0.17350。因此得到的图像为:

程序:

x1=[1.24,1.36,1.38,1.38,1.38,1.40,1.48,1.54,1.56];y1=[1.27,1.74,1.64,1.82,1.90,1.70,1.82,1.82,2.08];

x2=[1.14,1.18,1.20,1.26,1.28,1.30];y2=[1.78,1.96,1.86,2.00,2.00,1.96];

x3=1:0.01:2;

y3=1.4606*x3-0.1735/1.4356; plot(x1,y1,... '*',... x2,y2,... '+',... x3,y3);

legend('Af','Apf','分界线')

d3.m

x1=[1.24,1.36,1.38,1.38,1.38,1.40,1.48,1.54,1.56];y1=[1.27,1.74,1.64,1.82,1.90,1.70,1.82,1.82,2.08];

x2=[1.14,1.18,1.20,1.26,1.28,1.30];y2=[1.78,1.96,1.86,2.00,2.00,1.96];

zxec(x1,y1,1);zxec(x2,y2,1) x3=1:0.01:2;

y3=1.4606*x3-0.1735/1.4356; plot(x1,y1,... '*',... x2,y2,... '+',... x3,y3);

legend('Af类','APf类','分界线')

32.82.62.42.221.81.61.4 1AfApf分界线 1.11.21.31.41.51.61.71.81.92运行结果:

精度分析

分别将(1.24,1.80),(1.28,1.84),(1.40,2.04)代入方程

g(P(x,y))2.0968x1.4356y0.1735得:-0.1575,-0.1311,-0.1666

故这三个应属于Apf类的

6.4 (1) 利用

逼近多项式。

Legendre多项式构造函数f(x)x4在区间[0,1]上的二次最佳平方

*(x),(2) 设f(x)5x32x2x1,求一个次数不超过2的多项式p2使它成

为f(x)在区间[1,1]上的最佳一致逼近多项式.

(3) 利用截断Chebyshev级数求f(x)cosx在区间[0.5,1]上的近似三次最佳一致逼近多项式.

(4) 利用缩短幂级数法通过f(x)ex的5次Taylor多项式,求f(x)ex的三次近似多项式,并估计误差.

(1)调用Legendre函数

f=inline('x^4'); Legendre(f,2) ans =

0.171429*t^2 + 0.142857

(2)同(1)

f=inline('5*x^3+2*x^2-x+1'); Legendre(f,2) ans =

0.4*t^2 + 0.666667*t + 1.53333

(3)f=inline('cos(x)'); Chebyshev(f,3) (4)T=cat(1,6);

T=sym(T); syms t; T(1)=1; T(2)=t;

%求Chebyshev的各次多项式 for k=2:5

T(k+1)=collect(2*t*T(k)-T(k-1)); end

%taylor展开 format rat;

p3=taylor(exp(t),4), p4=taylor(exp(t),5); p5=taylor(exp(t),6); %展开多项式的对应系数 p5_xs=sym2poly(p5);

p4_xs=sym2poly(p4);

%展开多项式的对应的最高次项系数 a=p5_xs(1); b=p4_xs(1); %三次近似多项式 a*2^(-4)*T(6); b*2^(-3)*T(5);

p3_jsdxs=collect(p5-a*2^(-4)*T(6)-b*2^(-3)*T(5),x)

p3 =

t^3/6 + t^2/2 + t + 1 p3_jsdxs =

(17*t^3)/96 + (13*t^2)/24 + (383*t)/384 + 191/192

7.3设X为标准正态随机变量,即Xu0.1,0.2,0.3,N(0,1)。现分别取

,2.9,3.0,试设计算法计算30个不同的概率值PXua,

并将计算结果与概率论教科书中的标准正态分布函数表作比较。

(提示:PXua设n10

for a=0.1:0.1:3

b=4;n=10;h=(b-a)/n;

p=0;q=0;f=inline('exp((-x^2)/2)'); for i=0:n-1

p1=f((a+i*h+a+(i+1)*h)/2); p=p+p1; end

for j=1:n-1 q1=f(a+j*h); q=q+q1; end

format long

s=(1/sqrt(2*pi))*(h/6)*(f(a)+4*p+2*q+f(b)) end

12uaex22dx124uaex22dx)

结果: 求得的概率值与准确值的比较 标准正态分布函数表 复化Simpson公式计算结果 误差(10) 0.327005362270039 0.334323447380158 0.339255281809892 0.341958668619857 0.342716284110201 0.341892346639727 40.460172162722971 0.420740290560897 0.382088577811047 0.344578258389676 0.308537538725987 0.274253117750074 0.460139462186744 0.420706858216159 0.382054652282866 0.344544062522814 0.308503267097576 0.274218928515410 0.241963652223073 0.211855398583397 0.184060125346760 0.158655253931457 0.135666060946383 0.115069670221708 0.096800484585610 0.080756659233771 0.066807201268858 0.054799291699558 0.044565462758543 0.035930319112926 0.028716559816002 0.022750131948179 0.017864420562817 0.013903447513499 0.010724110021676 0.008197535924596 0.006209665325776 0.004661188023719 0.003466973803041 0.002555130330428 0.001865813300384 0.001349898031630

0.241929663311644 0.211821687826475 0.184026733797398 0.158622191452687 0.135633313093192 0.115037205509040 0.096768261388249 0.080724631768960 0.066775324333546 0.054767524041190 0.044533769129462 0.035898671181465 0.028684936179116 0.022718517515272 0.017832805567661 0.013871826372685 0.010692480205813 0.008165896953114 0.006178017945395 0.004629533589453 0.003435313851058 0.002523466314120 0.001834146452943 0.001318229317849 0.339889114290004 0.337107569219897 0.333915493619941 0.330624787699962 0.327478531910175 0.324647126680039 0.322231973610004 0.320274648109992 0.318769353120085 0.317676583680010 0.316936290810019 0.316479314609944 0.316236368860003 0.316144329070012 0.316149951559985 0.316211408140001 0.316298158630005 0.316389714819995 0.316473803810002 0.316544342659998 0.316599519830001 0.316640163080003 0.316668474410000 0.316687137810001

通过图像查看求得的概率值与准确值的比较:

x=0.1:0.1:3;

y=1-normcdf(x,0,1);

x1=[0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3]'; y1=[0.460139462186744 0.344544062522814 0.241929663311644 0.158622191452687 0.096768261388249 0.054767524041190 0.028684936179116 0.013871826372685 0.006178017945395 plot(x,y,'*r',... x1,y1,'.b')

0.420706858216159 0.308503267097576 0.211821687826475 0.135633313093192 0.080724631768960 0.044533769129462 0.022718517515272 0.010692480205813 0.004629533589453

0.382054652282866 0.274218928515410 0.184026733797398 0.115037205509040 0.066775324333546 0.035898671181465 0.017832805567661 0.008165896953114 0.003435313851058

0.002523466314120 0.001834146452943 0.001318229317849]';

legend('准确值','计算值')

0.50.450.40.350.30.250.20.150.10.050 00.511.522.53准确值计算值

从计算的积分值和准确值的散点图和误差,我们可以看出步长选为10步,运用辛卜生公式的精度非常高

7.5

xexdx, 分别用下列方法计算积分I0(1x)21 (1) Romberg算法;(2) 8点Gauss公式;

(3) 将积分区间分为5等份,用复化两点Gauss公式并比较结果。

(1) Romberg算法;

Roberg('f1',0,1) ans =

0.500000000000000

(2) 8点Gauss公式

查表得,8点gauss公式的节点和系数为:

x=[0.9602899 -0.9602899 0.7966665 -0.7966665 0.5255324 -0.5255324 0.1834346 -0.1834346 ];A=[0.102285 0.102285 0.2223810 0.2223810 0.3137067 0.3137067 0.3626838 0.3626838]; gauss8('f1',0,1,x,A)

7.6设某城市男子的身高X~N(170, 36)(单位:cm),应如何选择公共汽车门的

高度H使男子与车门碰头的机会小于1%。

(1)

%这里令n=10

for a=180:1:188

b=194;n=10;h=(b-a)/n;

p=0;q=0;f=inline('exp((-(x-170)^2)/72)'); for i=0:n-1

p1=f((a+i*h+a+(i+1)*h)/2); p=p+p1; end

for j=1:n-1 q1=f(a+j*h); q=q+q1; end

format long

P=(1/(6*sqrt(2*pi)))*(h/6)*(f(a)+4*p+2*q+f(b)) end

P = 0.047758637464544 P =0.033344869800652 P =0.022718517515272 P =0.015098521327612 P =0.009783695738731 P =0.006178017945395 P =0.003798722285391 P =0.002271601047586 P =0.001318229317849

(2)

for a=170:1:194 b=194;n=10;h=(b-a)/n;

p=0;q=0;f=inline('exp((-(x-170)^2)/72)'); for i=0:n-1

p1=f((a+i*h+a+(i+1)*h)/2); p=p+p1; end

for j=1:n-1 q1=f(a+j*h); q=q+q1; end

format long

P=(1/(6*sqrt(2*pi)))*(h/6)*(f(a)+4*p+2*q+f(b)); if P<=0.01 a break; end end a = 184

(3)

R=normrnd(170,6,10000,1);

j=0;

for i=1:10000 if R(i)>180 j=j+1; end end

disp('模拟数中大于180的数所占比例为:') j/10000

模拟数中大于180的数所占比例为: ans =

0.047600000000000

结果:模拟数中大于180的数所占比例为:

ans = 0.046300000000000

7.7

根据门特卡罗算法,我们在D所表示的圆形区域内任取n个随机数,然后分别计算

sinln(xy1)在每个随机数处的值,然后将这些数相加取其平均值。

程序及运行结果:

n=10000;

x=0:1;y=0:1;

fun=inline('sin(sqrt(log(x+y+1)))'); n0=0;

xh=x(2)-x(1); yh=y(2)-y(1); for i=1:n a=xh*rand+x(1); b=yh*rand+y(1);

f=feval(fun,a,b); if f<=0 n0=n0+1; end end

y=n0/n*xh*yh

输出结果为

y =0

8.2列出函数f(x)1lng(x)g(x)在区间[0, e]上的函数值表,并作出它的图

象;其中g(x)是初值问题

dgln[g(x)]2x,的解。 dxg(0)0.

建立m文件 mywork8_2

function dy=mywork8_2(x,y) dy=-2*x;

disp('求得的函数值为:')

[x y]=ode23(@mywork8_2,[0,exp(1)],0)plot(x,y); grid on;

title('f(x)的函数图象');

输出结果为:

求得的函数值为: x =

0 0.271828182845905 0.543656365691809 0.815484548537714 1.087312731383618 1.359140914229523 1.630969097075427 1.902797279921332 2.174625462767236 2.446453645613141 2.718281828459046 y =

0 -0.073890560989307 -0.295562243957226

-0.665015048903759 -1.182248975828904 -1.847264024732663 -2.660060195615034 -3.620637488476019 -4.728995903315616 -5.985135440133827 -7.389056098930652

函数图像为:

f(x)的函数图象0-1-2-3-4-5-6-7-800.511.522.53

8.8

鱼雷击舰问题。一敌舰在某海域内沿正北方向航行时,我方战舰恰位于敌舰的正西方1海里处。我舰向敌舰发射鱼雷,敌舰速度为0.42海里/分,鱼雷速度为敌舰速度的两部。试问敌舰航行多远时将被击中?

提示:本问题归结为常微分方程初值问题

1y'2,x(0,1),y''' 2(1x)y|x00,y'|x00.在计算机上用 MATLAB 求解 解鱼雷运动轨迹的参数方程 因为假设 v0 =1,所以敌舰所处位置为(1,t) 。

建立 m-文件 eq2.m 如下: function dy=eq2(t,y) dy=zeros(2,1);

dy(1)=2*(1-y(1))/sqrt((1-y(1))^2+(t-y(2))^2); dy(2)=2*(t-y(2))/sqrt((1-y(1))^2+(t-y(2))^2); 取 t0=0,tf=2,建立主程序

[t,y]=ode45('eq2',[0 2],[0 0]); Y=0:0.01:2; plot(1,Y,'-')

hold on, plot(y(:,1),y(:,2),'*')

2.521.510.5000.20.40.60.811.21.4

如图,在 x=1,y≈0.667 左右相交。

接着,我们运用二分法修改 tf 的值

[t,y]=ode45('eq2',[0 0.667],[0 0]); Y=0:0.01:2; plot(1,Y,'-') hold on, plot(y(:,1),y(:,2),'*')

21.81.61.41.210.80.60.40.2000.20.40.60.811.21.4

即鱼雷与敌舰相遇地点为(1,0.667) ,所以相遇所耗时间 t≈88.9 秒

附录:作业中所使用函数

因篇幅问题不能全部显示,请点此查看更多更全内容