2023年数值计算大作业

济树帆力挺）、孝

课程设计

课程名称：___________________________________

设计题目：___________________________________

学号：___________________________________

姓名:___________________________________

完毕时间:

题目一：非线性方程求根

一摘要

非线性方程的解析解通常很难给出，因此非线性方程的数值解就尤为重要。本实

验通过使用常用的求解方法二分法和Newton法及改善的Newton法解决几个

题目,分析并总结不同方法解决问题的优缺陷。观测迭代次数，收敛速度及初值选

取对迭代的影响。

用Newton法计算下列方程

⑴x3-x-l=0,初值分别为%=1,%=0.45,%=0.65；

（2）丁+94/—389x+294=0其三个根分别为1,3,-98。当选择初值%=2时

给出结果并分析现象，当£=5x10-6,迭代停止。

解：1）采用MATLAB进行计算；

一方面定义了Newton法：

functionkk=newton（fzdf,xOzto1,N）

%NewtonMethod（牛顿法）

%Thefirstparameterfisaexternalfunctionwithrespectt

oviablex.（第一个参数也就是本题所用的函数f）

%Thesecondparameterdfisthefirstorderdiffentialfunc

tionoffx.（第二个参数也就是本体所用函数f的导数方程df）

%x0isinitialiterationpoint（初值）.

%tolisthetoleranceofthe1oop（精度）.

%Nisthemaximumnumberofiterations（循环上限）.

x=x0；

f0=eval（f）；df0=eva1（df）；

n=0;

disp('[nxnxn+1fn+1],);

whilen<=N

xl=xO-fO/df0;

x=x1;

f1=eva1(f);

X=[nzx0,xlzf1];

disp(X);

ifabs(x0-xl)<to1

fprintf(*Theprocedurewassuccessful.*)

kk=X;

return

else

n=n+l;

x0=xl；f0=f1;

end

end

ifn==N+l

fprintf(*themethodfailedafterNiterations.

kk=0;

End

我们把Newton法存为.m格式的文献；

之后我们运营程序：

clear；

clc;

symsx

f=xA3-x-l;

df=diff(f,x);

x=newton(f,df,1,0.0001,50);

x

会得到一下结果

Enxnxn+1fn+1]

01.00001.50000.8750

1.00001.50001.0625-0.8630

2.00001.06251.49400.8408

到第50次迭代时候会出现该问题:

47.00001.48981.0814-0.8167

48.00001.08141.48980.8167

49.00001.48981.0814-0.8167

50.00001.08141.48980.8167

themethodfailedafterNiterations.

x=

0；

同样测试x0=0.45、0.65得不出结果，判断出初值离真值太远，所以我们采用

牛顿下山法进行计算迭代：

我们定义了其中的f函数和df函数，并且分别存为.m格式的文献,其代码如下:

f:

functiony=f(x)

y=xA3-x-1;

df:

functiony=df(x)

y=3*xA2-l；

之后我们定义newton下山法同时也存为.m的程序：

function[x,i]=downnewton(fzdfzxO,tol)

k=0;

i=l;

disp(*[nxnxn+1fn+1]1);

while(k==0)

fx=feval('f/,xO)；

7

dfx=feval(*dfzxO);

t=0;

u=1；

while(t==0)

dx=-fx/dfx；

xl=x0+u*dx；

fxl=feva1('f*zx1)；

fx0=feval(1f',xO)；

if(abs(fx1)>abs(fxO));

u=u/2;

e1se

t=l;

end

end

X=[i,xO,xl,fxl];

disp(X);

if(abs(fxl)<tol)

k=l;

eIse

x0=x1;

i=i+1;

end

end

x=xl;

i=i；

end

之后带入X0=0.45；

downnewton(T,*df,0.45,10A(-6))

[nxnxn+1fn+1

1.00000.4500-0.4155-0.6562

2.0000-0.4155-0.5857-0.6152

3.0000-0.5857-0.5754-0.6151

4.0000-0.5754-0.5782-0.6151

5.0000-0.5782-0.5773-0.6151

6.0000-0.5773-0.5774-0.6151

7.0000-0.5774-0.5773-0.6151

8.0000-0.5773-0.5774—0.6151

9.0000-0.5774-0.5774-0.6151

10.0000-0.5774-0.5774-0.6151

11.0000-0.57741.3131-0.0490

12.00001.31311.32480.0005

13.00001.32481.32470.0000

ans=

1.3247

带入x0=0.6;

downnewton('f,zdf,0.6,10A(-6))

[nxnxn+1fn+1]

1.00000.60001.1406-0.6566

2.00001.14061.36680.1866

3.00001.36681.32630.0067

4.00001.32631.32470.0000

5.00001.32471.32470.0000

ans

1.3247

带入xO=l;

downnewton('f,*df,1,10A(-6))

[nxnxn+1fn+l]

1.00001.00001.50000.8750

2.00001.50001.34780.1007

3.00001.34781.32520.0021

4.00001.32521.32470.0000

ans=

1.3247

2)同样采用Newton下山法：

重新定义f、df:

f:functiony=f(x)

y=xA3+94*xA2-389*x+294;

df：functiony=df(x)

y=3*xA2+188*x-389；

再带入初值x0=2；

downnewton,2,5*10"(-6))

nxnxn+1fn+1

12-980

ans

-98

得出x=-98；

分析：先画出该函数的图像；

x=(-100:.1:100);ezplot('xA3+94*x2-389*x+294100100J)

得出该图像如图：

根据牛顿法的几何解释，在x0=2的点做切线，与y相交，交点的横坐标值为x=-98则结束了

该现象。

题目二:线性方程组求解

一摘要

对于实际的工程问题，很多问题归结为线性方程组的求解。本实验通过实际题目

掌握求解线性方程组的数值解法，直接法或间接法。

有一平面机构如图所示，该机构共有13条梁（图中标号的线段）由8个较接点（图

中标号的圈）联结在一起。上述结构的1号钱接点完全固定，8号较接点竖立方

向固定，并在2号、5号和6号钱接点，分别有如图所示的10吨、15吨和20

吨的负载，在静平衡的条件下，任何一个较接点上水平和竖立方向受力都是平衡

的，以此计算每个梁的受力情况。

»1015*20

令。=1/V2,假设f为各个梁上的受力，例如对8号钱接点有

九+舫2=。

对5号钱接点，则有

05+4=班+/|0

次+6+次=15

针对各个钱接点，列出方程并求出各个梁上的受力。

解:针对此题我们采用雅克比迭代法；

一方面我们先写出Jacobi迭代的程序，并且存为.m的形式：

function[x,n]=jacobi(A,b,x0,eps,varargin)

ifnargin==3

eps=1.Oe-6；

M=200；

elseifnargin<3

error

return

elseifnargin==5

M=varargin{1}；

end

D=diag(diag(A));

L=-tril(A,-1);

U=-triu(Ar1)；

B=D\(L+U);

f=D\b;

x=B*xO+f；

n=l;

whilenorm(x-xO)>=eps

xO=x;

x=B*xO+f;

n=n+1；

if(n>=M)

disp（'Warning:迭代次数太多，也许不收敛!,）;

return;

end

end

之后我们根据节点进行计算杆的力，设受拉为正；

1）由于角度为45。,所以正弦值和余弦值相等都设为a=2N-1/2）;

则可列方程:afi=O;

afi+f2=0；

f3=10;

f6=0;

afi+f3+af5=O；

afi-f4-af5=0;

f4-f8=0；

f7=o；

af5+f7+af9=15;

储=20；

flO-f13=0；

afi2=0;

afi2+fi3=0;

输入到mat1ab中有如下:

A=[2A(-1/2)000000000000;2-(・l/2)100000000000;00

10000000000;01000-10000000;2"(-l/2)0102人(-1/2)0

0000000;2^(-l/2)00-1一2八(一1/2)00000000;00

01000-100000;0000001000000；00002A(-1/2)01

02A(-l/2)0000:0000000000100；000000000

100-l;000000000002A(-1/2)0;000000000002人(-1/2)

1;],b=[001000000152000031

得出矩阵A和b

A=

0.707100。0。00。000OoO0