数值分析 - 编程作业(4)

y0(i)=1/(8.973276262446161e-003+8.416900197047816e-004*1/x0(i)); end;

plot(x0,y0,x,y); 如下图：

第七章 数值微分与数值积分

1、考纽螺线的形状像钟表的发条，也称回旋曲线，它在直角坐标系中的参数方程为

s12?x(s)?cos?02atdt ??s1?y(s)??sinat2dt02?曲线关于原点对称。取a=1，参数s的变化范围[?5,5]，容许误差限分别是10?3和

10?7。选择适当的节点个数，利用数值积分方法计算曲线上点的坐标，并画出曲

线的图形。

解：由于题目中规定了误差限，又要考虑到复杂度及运行时间，故选择误差限为

10?3。可以选用梯形公式的逐次分半算法，并且s的点数不宜取过多，本例中一

共取101个点。 程序如下：

先建两个M文件，保存函数x与y。M文件如下： 文件xf.m保存x?s? function f1=xf(t)

f1=cos(t^2/2); 文件yf.m保存y?s?

function f2=yf(t) f2=sin(t^2/2);

用复化梯形公式计算积分，程序如下： function [x,y]=fuhuatixing() dx=10^-3;dy=10^-7; i=1;

for s=-5:0.1:5 m=1; b=s; a=0;

h=(b-a)/2;

T=h*(xf(a)+xf(b)); T0=5;

while abs(T-T0)>3*dx

Fx=0; T0=T; for k=1:2^(m-1)

Fx=Fx+xf(a+(2*k-1)*h); end

T=T0/2+h*Fx; m=m+1; h=h/2; end x(i)=T; i=i+1; end j=1;

for s=-5:0.1:5 n=1; b=s; a=0;

h=(b-a)/2;

T=h*(yf(a)+yf(b)); T0=5;

while abs(T-T0)>3*dy Fy=0; T0=T; for k=1:2^(n-1)

Fy=Fy+yf(a+(2*k-1)*h); end

T=T0/2+h*Fy; n=n+1; h=h/2; end y(j)=T; j=j+1; end

plot(x,y,'o-g');

在命令窗口中输入[x,y]=fuhuatixing()

可得x与y的值，即曲线上点的坐标，在本例中s从-5到5步长为0.1计算了101个点的坐标值，得到曲线图形如下图所示：