数值分析大作业三四五六七完整版.docx

资源描述

数值分析大作业三四五六七完整版.docx

《数值分析大作业三四五六七完整版.docx》由会员分享，可在线阅读，更多相关《数值分析大作业三四五六七完整版.docx（25页珍藏版）》请在冰豆网上搜索。

数值分析大作业三四五六七完整版.docx

数值分析大作业三四五六七完整版

数值分析大作业三四五

六七

HENsystemofficeroom[HEN16H-HENS2AHENS8Q8-HENH1688]

大作业三

I.给定初值Xo及容许误差£,编制牛顿法解方程f{x）=0的通用程序.

解：

Matlab程序如下：

函数m文件：

functionFu=fu（x）

Fu=x"3/3-x;

end

函数m文件：

functionFu=dfu（x）

Fu=x*2-1;

end

用Newton法求根的通用程序

clear;

x0=inputC请输入初值x0:

‘）；

ep=input（，请输入容许i吴崔：

’）；

flag=l;

whileflag==l

xl=x0-fu（xO）/dfu（xO）;

ifabs（xl-xO）

flag=0;

end

x0=xl;

end

fprintf（'方程的•个近似解为:

%f\n',xO）;

寻找最大6值的程序：

clear

eps=input（'请输入搜索淸度:

;

ep=input（，诘输入容许谋差：

’）;

flag=l;

k=0;

x0=0;

whileflag==l

sigma=k*eps;

x0=sigma;

k=k+l;

m=0;

flagl=l;

whileflagl==l&&m〈=10"3

xl=xO-fu（xO）/dfu（xO）;

ifabs（xl-xO）

flagl=O;

end

m=m+l;

xO=xl;

end

ifflagl==labs（xO）>=ep

flag=O;

end

fprintf赧人的sigma值为：

%f\n*,sigma）;

求下列方程的非零根

解：

Matlab程序为:

（1）主程序

clear

clcformatlongx0=765;

N=100;

errorlim=10"（-5）;

x=xO-f（xO）/subs（df（）,xO）;n=l;

whilen

x=xO-f（xO）/subs（df0,xO）;ifabs（x-xO）>errorlimn=n+l;

else

break;

end

xO=x;

end

disp（[‘迭代次数：

n=*,num2str（n）J）

disp（['所求IE零根:

正根xl=',num2str（x）,'负根x2=',num2str（-x）J）

（2）子函数非线性函数f

functiony=f（x）

y=log（（513+*x）/*x））-x/（1400*;

end

（3）•函数非线性函数的•阶导数df

functiony=df（）

symsxl

y=log（（513+*xl）/*xl））-xl/（1400*;

y=diff（y）;

end

运行结果如下：

迭代次数:

n=5

所求非零根：

正根xU负根x2=

大作业四

试编写MATLAB函数实现Newton插值，要求能输

出插值多项式.对函数/*（%）=―在区间[-5,5]1+4Q

上实现10次多项式插值.

分析：

（1）输出插值多项式。

（2）在区间[-5,5]内均匀插入99个节点，计算这些节点上函数f（x）的近似值，并在同一张图上画出原函数和插值多项式的图形。

（3）观察龙格现象，计算插值函数在各节点处的误差，并画出误差图。

解：

Matlab程序代码如下：

%此函数实现y=l/（l+4*x*2）的n次Newton插值，n由调用函数时指定

%函数输出为插值结果的系数向量（行向量〉和插值多项式

function[ty]=func5（n）

xO=linspace（-5,5,n十1）';

y0=l./（l.+4.*x0."2）;

b=zeros（1,n+1）;

fori=l:

n+l

s=0;

forj=l:

t=l；

fork=l:

ifk~=j

t=（xO（j）-xO（k））*t;

end;

s=s+yO（j）/t;

end;

b（i）=s;

end;

t=linspace（0,0,n+1）;

fori=l:

s=linspace（0,0,n+1）;

s（n+l-i:

n+1）=b（i+1）.*poly（x0（l:

i））;

t=t+s;

end;

t（n+l）=t（n+l）+b（l）;

y=poly2sym（t）;

10次插值运行结果：

[bY]=func5（10）

b二

Columns1through4

Columns5through8

Columns9through11

Y二

b为插值多项式系数向量，Y为插值多项式。

插值近似值：

xl=linspace（-5,5,101）;

x=xl（2:

100）;

y=polyval（b,x）

y二

Columns1through12

绘制原函数和拟合多项式的图形代码:

plot（x,1./（1+4.*x・"2））

holdall

plot（x,y,'r'）

xlabel（'X’）ylabel（'Y'）title（'Runge现象'）gtextC原函数'）gtext十次牛顿插值多项式'）詁差计数并绘制误差图：

holdoff

ey=l./（1+4・*x."2）-y

ey=

Columns1through12

Columns13through24

Columns25through36

Columns37through48

Columns49through60

Columns61through72

Columns73through84

Columns85through96

plot（x,ey）

xlabelfX'）

ylabel（'ey'）

title（'Runge现象误差图'）

输出结果为：

大作业五

解：

Mat1ab程序为：

x=[-520,-280,,-7&,,0,,,78,,280,520]';

y=[0,-30,-36,-35,,,0,,,35,36,30,0]';

n=13;

喘求解M

fori=1:

n-l

h（i）=x（i+l）-x（i）;

end

fori=2:

n-l

a（i）=h（i-l）/（h（i-l）+h（i））;

b（i）=l-a（i）;

c（i）=6*（（y（i+l）-y（i））/h（i）-（y（i）-y（i-l））/h（i-l））/（h（i-l）+h（i））;

end

a（n）=h（n-l）/（h（l）+h（n-l））;

b（n）=h（l）/（h（l）+h（n-l））;

c（n）=6/（h（l）+h（n-l））*（（y⑵-y（l））/h⑴-（y（n）-y（n-l））/h（n-l））;

A（l,1）=2;

A（l,2）=b

（2）;

A（l,n-1）=a

（2）;

A（n-1,n-2）=a（n）;

A（n-1,n-1）=2;

A（n-1,1）=b（n）;

fori=2:

n-2

A（iti）=2;

A（i,i+1）=b（i+l）;

A（i,i-1）=a（i+l）;end

C=c（2:

n）;

m=A\C;

M（l）=m（n-l）;

M（2:

n）=m;

xx=-520:

520;

fori=1:

forj=1:

n-l

ifx（j）<=xx（i）&&xx（i）

break;

end

yy（i）=M（j+l）*（xx（i）-x（j）厂3/（6*h（j））-M（j）*（xx（i）-x（j+1））*3/（6*h（j））+（y（j+l）-M（j+1）*h（j）*2/6）*（xx（i）-x（j））/h（j）-（y（j）-M（j）*h（j）"2/6）*（xx（i）-x（j+l））/h（j）;end;

fori=52:

101

yy（i）=-yy（102-i）;

end;

fori=1:

xx（i）=-xx（i）;

end

plot（xx,yy）;

holdon;

fori=1:

n/2

x（i）=-x（i）;

end

plot（x,y,'bd‘）；

titleC机翼外形曲线'）;

输出结果：

运行文件，得到

（1）编制求第一型3次样条插值函数的通用程序；

（2）已知汽车门曲线型值点的数据如下：

解：

（1）Matlab编制求第一型3次样条插值函数的通用程序：

function[Sx]=Threch（X,Y,dyO,dyn）

%X为输入变量x的数值

喘Y为函数值y的数值

瓯dyO为左端•阶导数值

%dyn为右端"阶导数值

瓯Sx为输出的函数农达式

n=length（X）~l;

d=zeros（n+1,1）;

h=zeros（1,n-1）;

fl=zeros（1,n-1）;

f2=zeros（l,n-2）;

fori=l:

n喘求函数的•阶筮瀚

h（i）=X（i+1）-X（i）;

fl（i）=（Y（i+l）-Y（i））/h（i）;

end

fori=2:

n嗚求函数的二阶签商

f2（i）=（f1（i）-fl（i-1））/（X（i+1）-X（i-1））;

d（i）=6*f2（i）;

end

d（l）=6*（fl（l）-dy0）/h（l）;

d（n+l）=6*（dyn-fl（n-l））/h（n~l）;

%赋初值

A=zeros（n+1,n+1）;

B=zeros（1,n-1）;

C=zeros（1,n-1）;

fori=l:

n-l

B（i）=h（i）/（h（i）+h（i+l））;

C（i）=l-B（i）;

end

A（l,2）=1；

A（n+l,n）=1;

fori=l:

n+l

A（i,i）=2;

end

fori=2:

A（i,i-1）=B（i-l）;

A（ifi+1）=C（i-l）;

end

M=A\d;

symsx;

fori=l:

Sx（i）=collect（Y（i）+（fl（i）-（M（i）/3+M（i+l）/6）*h（i））*（x-X（i））+M（i）/2*（x-X（i））*2+（M（i+l）-M（i））/（6*h（i））*（x-X（i））*3）;

digits（4）;

Sx（i）=vpa（Sx（i））;

end

fori=l:

dispCS（x）='）;

fprintf（'%s（%d,%d）\n,char（Sx（i））,X（i）,X（i+1））;

end

S=zeros（1,n）;

fori=l:

x=X⑴十；

S（i）=Y（i）+（f1（i）-（M（i）/3+M（i+1）/6）*h（i））*（x-X（i））+M（i）/2*（x-X（i））*2+（M（i+1）-M（i））/（6*h（i））*（x-X（i）厂3;

end

disp（*S（i+‘）;

dispCiX（i+S（i+'）;

fori=l:

fprintfC%d%.4f%.4f\n,i,X（i）+,S（i））;

end

输出结果：

»X=[012345678910];

»Y=[];

»Threch（X,Y,,

S（x）=

—

*x"2

-*x"3+

（0,1）

S（x）=

—

*x"2

-*x"3+

（1,2）

S（x）二

—

*x"2

-*x"3+

（2,3）

S（x）=

一*x

-*x"3+

（3,4）

S（x）=

—

*x"2

+*x"3-

（4,5）

S（x）二

一*x

-*x"3+

（5,6）

S（x）=

—

*x"2

+*x"3-

（6,7）

S（x）=

一*x

-*x"3+

（7,8）

S（x）=

—

*x"2

+*x"3-

（&9）

S（x）二

—

*x"2

+*x"3-

（9,10）

S（i+

•

X（i+

S（i+

10ans

[-*x"3-*x"2+*x+，-*x"3-*x"2+*x+,-*x"3-*x"2+*x+,-*x"3+*x"2-*x+,*x"3-*x"2+*x-,-*x"3+*x"2-*x+,*x"3-*x"2+*x-,-*x"3+*x"2-*x+,*x"3-*x"2+*x-,*x"3-*x"2+*x一]

大作业六

1、炼钢厂出钢时所用的圣刚睡的钢包，在使用过程中111于钢液及炉渣对包衬耐火材料的侵蚀，使其容积不断增大，经试验，钢包的容积与相应的使用次数的数据如下：

（使用次数X,容积y）

f/?

*1/x

选用

双曲

线

l/y=a

对使用最

乘法

数据

进行

拟合。

解：

Matlab程序如下：

functiona=nihehanshu（）

x0=[2356791011121416171920];

y0=[1;

A=zeros（2,2）;

B=zeros（2,1）;

a=zeros（2,1）;

x=l・/x0;

y=l./y0;

A（l,1）=14;

A（l,2）=sum（x）;

A（2,1）=A（1,2）;

A（2,2）=sum（x.*2）;

B（l）=sum（y）;

（2）=sum（x.*y）;

a=A\B;

y=l./（a（l）+a

（2）*1./x0）;subplot（1,2,2）;

plot（xO,yO-y,'bd-'）;titleC拟合曲线误差'）;subplot（1,2,1）;

plot（xO,yO,'go'）;

holdon;

x=2:

20;

y=l./（a（l）+a

（2）*1./x）;plot（x,y,'r*-'）;

legend（'散点,，'拟介曲线图1/y=a

（1）+a

（2）*l/x'）;title（'最小二乘法拟介曲线'）;

试验所求的系数为：

nihehanshu

ans=

则拟合曲线为丄=0.008973276262446+0.000841690019705-

拟合曲线图、散点图、误差图如下：

2、下面给出的是乌鲁木齐最近1个月早晨7:

00左右（新媼时间）的天气预报所得到的温度，按照数据找出任意次曲线拟合方程和它的图像。

用MATLAB编程对上述数据进行最小二乘拟合。

2008年10月--11月26日

温度

天数

温度

盍

温度

解：

Matlab的程序如下：

x=[l:

30];

y=[9,10,11,12,13,14,13,12,11,9,10,11,12,13,14,12,11,10,9,8,7,&9,11,9,7,6,5,3,1];

al=polyfit（x,5%3）%三次多项式拟合％

a2=polyfit（x,y,9）%九次多项式拟合%

a3=polyfit（x,y,15）%十五次多项式拟合%

bl=polyval（al,x）

b2=polyval（a2,x）

b3=polyval（a3,x）

rl=sum（（y-bl）.'2）滋三次多项式i吴差平方和％

r2=sum（（y-b2）.*2）%九次次多项式误差平方和％

r3=sum（（y-b3）."2）弔十五次多项式误差平方和％

plot（x,y/**）%用*画出x,y图像%

holdon

plot（x,bl,'r）%用红色线画出x,bl图像%

holdon

plot（x,b2,'g）%用绿色线画出x,b2图像%

holdon

plot（x,b3,'b:

o‘）弘用蓝色o线画出x,b3图像%

试验结果为：

al=

Columns1through2

Columnsa2=

Columns

Columnsa3=

Columns

Columnsbl=

Columns

Columnsb2=

Columns

through4

through2

through4through6through8through10

through2

through4

through6

through8through10

through12

through14

through16

through2

through4

through6

through8through10

through12

through14

through16

through18

through20

through22

through24

through26

through28

through30

through2

through4

through6

through8through10

through12

through14

through16

through18

through20

through22

through24

Columns

throughthroughthrough

b3=

Columns1through2

Columns

through4

Columns

through6

Columns

through8

Columns

through10

Columns

through

Columns

through

Columns

through

Columns

through

Columns

through

Columns

through

Columns

through

Columns

through

Columns

through

Columns

through

rl=r2=r3=

其中的最后图像为:

大作业七

用Romberg求积法计算积分

的近似值，要求误差不超过①5x1旷

解：

Matlab程序为：

玄被积函数m文件：

functionFx=fx（x）

Fx=l/（l+100*x*x）;

end

%Romberg求积法计算积分的通用程序

functionRomberg0

clear;

a=input（'请输入积分下限:

a='）;b=input（'请输入积分I•.限：

b='）;eps=input（，请输入允许持度：

eps='）;

%=====计算Tn=====%

functionTn=T（n）

Tn=0;

h=（b-a）/n;

x=zeros（1,n+1）;

fork=l:

n+l

x（k）=a+（k-l）*h;

end

forj=l:

Tn=Tn+h*（fx（x（j））+fx（x（j+1）））/2;

end

%=======计算Sn======%

functionSn=S（n）

Sn=4/3*T（2*n）-l/3*T（n）;

end

%========计算Cn=====%

functionCn=C（n）

Cn=16/15*S（2*n）-1/15*S（n）;

end

%=====计算Rn=====%

functionRn=R（n）

Rn=64/63*C（2*n）-1/63*C（n）;

end

賀====d十算满足允许精度的Rn，并打印输｛!

«,========%

i=l；

flag=l；

whileflag==l

ifabs（R（2"i）-R（2"（i-l）））/255

flag=0

展开阅读全文