牛顿柯特斯 matlab,数值分析上机实验报告

数值分析上机实验报告"/>

牛顿柯特斯 matlab,数值分析上机实验报告

Y=zeros(1,n+1); Y(1)=0; for i=2:n+1 x=X(i-1); y=Y(i-1);

Y(i)=Y(i-1)+eval(F)*h; end Y(6) ans =

0.1186 Euler校正法

Y1=zeros(1,n+1); Y1(1)=0;

for i=2:n+1 x=X(i-1); y=Y1(i-1);

ty=Y1(i-1)+eval(F)*h; Y1(i)=Y1(i-1)+h/2*eval(F); x=X(i); y=ty;

Y1(i)=Y1(i)+h/2*eval(F); end Y1(6) ans =

0.1450 准确值 temp=[];

f=dsolve('Dy= x*x+x-y ','y(0)=1','x'); df=zeros(1,n+1); for i=1:n+1

temp=subs(f,'x',X(i)); df(i)=double(vpa(temp)); end

disp(' 步长 Euler法 准确值');

disp([X',Y',Y1',df']); 画图观察效果

运行上述程序得到结果如下:

步长 Euler法 Euler预测-校正公式 准确值 0 0 0 1.0000 0.1000 0 0.0055 0.9100 0.2000 0.0110 0.0219 0.8400 0.3000 0.0339 0.0501 0.7900 0.4000 0.0695 0.0909 0.7600 0.5000 0.1186 0.1450 0.7500 0.6000 0.1817 0.2130 0.7600 0.7000 0.2595 0.2954 0.7900 0.8000 0.3526 0.3929 0.8400 0.9000 0.4613 0.5059 0.9100

1.0000 0.5862 0.6348 1.0000 同时作出的对比图如下:

Euler预测-校正公式 Euler法和Euler预测-校正法解常微分方程10.90.80.70.60.50.40.30.20.10 00.10.20.30.40.50.60.70.80.91准确值Euler法Euler预测-校正法

由图看出Euler和Euler预测-校正法,都不好 2.利用Euler方法计算积分

?x0etdt在点x=0.5，1,1.5，2的近似值。

21四阶经典R-K公式作数值计算 clc;

F='3*y/(1+x)'; a=0; b=1; h=0.2; n=(b-a)/h; X=a:h:b;

Y=zeros(1,n+1); Y(1)=1; for i=1:n x=X(i); y=Y(i);

K1=h*eval(F); x=x+h/2; y=y+K1/2; K2=h*eval(F); x=x;

y=Y(i)+K2/2; K3=h*eval(F); x=X(i)+h; y=Y(i)+K3; K4=h*eval(F);

Y(i+1)=Y(i)+(K1+2*K2+2*K3+K4)/6; end 准确解 temp=[];

f=dsolve('Dy=3*y/(1+x)','y(0)=1','x'); df=zeros(1,n+1); for i=1:n+1

temp=subs(f,'x',X(i));

df(i)=double(vpa(temp)); end

disp(' 步长 四阶经典R-K法 准确值'); disp([X',Y',df']); %画图观察效果 figure;

plot(X,df,'k*',X,Y,'--r'); grid on;

title('四阶经典R-K法解常微分方程'); legend('准确值','四阶经典R-K法');

运行上述程序得到如下结果

步长 四阶经典R-K法 准确值 0 1.0000 1.0000 0.2000 1.2428 1.2428 0.4000 1.5836 1.5836 0.6000 2.0442 2.0442 0.8000 2.6510 2.6511

1.0000 3.4365 3.4366 作出的函数图形如下:

四阶经典R-K法解常微分方程3.5准确值四阶经典R-K法3 2.521.51 00.10.20.30.40.50.60.70.80.91

a=temp; b=b; end end;

disp('k a(k) b(k) x(k) f(x) '); H=[[1:k]',A',B',X',F']; disp(H); disp('x='); disp(X(k)); disp('y='); disp(yt);

format short

k a(k) b(k) x(k) 1.000000000000000 0.760287230697899

2.000000000000000 0.760287230697899

3.000000000000000 0.760287230697899

4.000000000000000 0.760287230697899

5.000000000000000 0.760287230697899

6.000000000000000 0.760287230697899

7.000000000000000 0.760287230697899

8.000000000000000 0.760287230697899

9.000000000000000 0.760287230697899

10.000000000000000 0.760287230697899

11.000000000000000 0.760287230697899

12.000000000000000 0.760287230697899

13.000000000000000 0.760287230697899

14.000000000000000 0.760287230697899

15.000000000000000 0.760287230697899

16.000000000000000 f(x)

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.000000000000000 1.000000000000000 0.500000000000000 0.500000000000000 0.500000000000000 0.500000000000000 0.500000000000000 0.500000000000000 0.500000000000000 0.500000000000000 0.500000000000000 0.500000000000000 0.500000000000000 0.500000000000000 0.500000000000000 0.500000000000000 0.500000000000000 0.500000000000000 1.000000000000000 1.000000000000000 1.000000000000000 1.000000000000000 1.000000000000000 1.000000000000000 1.000000000000000 1.000000000000000 1.000000000000000 1.000000000000000 1.000000000000000 1.000000000000000 1.000000000000000 1.000000000000000

0.760287230697899

17.000000000000000 0 1.000000000000000 0.500000000000000 0.760287230697899

18.000000000000000 0 1.000000000000000 0.500000000000000 0.760287230697899

19.000000000000000 0 1.000000000000000 0.500000000000000 0.760287230697899

20.000000000000000 0 1.000000000000000 0.500000000000000 0.760287230697899

21.000000000000000 0.760287230697899

22.000000000000000 0.760287230697899

23.000000000000000 0.760287230697899

24.000000000000000 0.760287230697899

25.000000000000000 0.760287230697899

26.000000000000000 0.760287230697899

27.000000000000000 0.760287230697899

28.000000000000000 0.760287230697899

29.000000000000000 0.760287230697899

30.000000000000000 0.760287230697899

31.000000000000000 0.760287230697899

32.000000000000000 0.760287230697899

33.000000000000000 0.760287230697899

34.000000000000000 0.760287230697899

35.000000000000000 0.760287230697899

36.000000000000000 0.760287230697899

37.000000000000000 0.760287230697899

38.000000000000000 0 1.000000000000000 0 1.000000000000000 0 1.000000000000000 0 1.000000000000000 0 1.000000000000000 0 1.000000000000000 0 1.000000000000000 0 1.000000000000000 0 1.000000000000000 0 1.000000000000000 0 1.000000000000000 0 1.000000000000000 0 1.000000000000000 0 1.000000000000000 0 1.000000000000000 0 1.000000000000000 0 1.000000000000000 0 1.000000000000000 0.500000000000000 0.500000000000000 0.500000000000000 0.500000000000000 0.500000000000000 0.500000000000000 0.500000000000000 0.500000000000000 0.500000000000000 0.500000000000000 0.500000000000000 0.500000000000000 0.500000000000000 0.500000000000000 0.500000000000000 0.500000000000000 0.500000000000000 0.500000000000000

0.760287230697899

39.000000000000000 0 1.000000000000000 0.500000000000000 0.760287230697899

40.000000000000000 0 1.000000000000000 0.500000000000000 0.760287230697899

41.000000000000000 0 1.000000000000000 0.500000000000000 0.760287230697899

42.000000000000000

实验报告六

一． 实验名称

方程求根

二． 实验目的

掌握二分法与Newton法的基本原理和应用。

三． 实验原理

对于二分法，其数学实质就是说对于给定的待求解的方程f(x)，其在[a,b]上连续，f(a)f(b)<0，且f(x)在[a,b]内仅有一个实根x，取区间中点c，若，则c恰为其根，否则根据f(a)f(c)<0是

*否成立判断根在区间[a,c]和[c,b]中的哪一个，从而得出新区间，仍称为[a,b]。重复运行计算，直至满足精度为止。这就是二分法的计算思想。

Newton法通常预先要给出一个猜测初值x0，然后根据其迭代公式

xk?1?xk?*f(xk) 'f(xk)*产生逼近解x的迭代数列{xk}，这就是Newton法的思想。当x0接近x时收敛很快，但是当x0选择不好时，可能会发散，因此初值的选取很重要。另外，若将该迭代公式改进为

xk?1?xk?rf(xk) 'f(xk)其中r为要求的方程的根的重数，这就是改进的Newton法，当求解已知重数的方程的根时，在同种条件下其收敛速度要比Newton法快的多。 程序设计：

本实验采用Matlab的M文件编写。其中待求解的方程写成function的方式，如下 function y=f(x); y=-x*x-sin(x);

写成如上形式即可，下面给出主程序。 二分法源程序： clear

b=1.5; a=0; R=1; k=0;

while (R>5e-6) ; c=(a+b)/2;

if f12(a)*f12(c)>0; a=c; else b=c; end R=b-a; k=k+1; end x=c

Newton法及改进的Newton法源程序： Clear;

f=input('请输入需要求解函数>>','s') df=diff(f); miu=2; %%%初始值x0

x0=input('input initial value x0>>');

k=0;%迭代次数 max=100;%最大迭代次数

R=eval(subs(f,'x0','x'));%求解f(x0)，以确定初值x0时否就是解 while (abs(R)>1e-8)

x1=x0-miu*eval(subs(f,'x0','x'))/eval(subs(df,'x0','x')); R=x1-x0; x0=x1; k=k+1;

if (eval(subs(f,'x0','x'))<1e-10); break end

if k>max;%如果迭代次数大于给定值，认为迭代不收敛，重新输入初值 ss=input('maybe result is error,choose a new x0,y/n?>>','s'); if strcmp(ss,'y')

x0=input('input initial value x0>>'); k=0;

else break end end end

k;%给出迭代次数

x=x0；

四． 实验内容

1． 用比例求根法求f(x)=1-xsinx=0在区间(0，1)内的一个根，直到近似根xk满足精度|f(xk)|<0.00005时终止计算。 解: clear format long x=1;

f=inline('1-x*sin(x)'); a=1; b=pi/2;

while abs(a-b)>0.00001 if f((a+b)/2)<0 b=(a+b)/2; x=b; end

if f((a+b)/2)>0 a=(a+b)/2; x=a; end

if f((a+b)/2)==0 x=(a+b)/2; break end end x

运行结果： x =

1.114153168596455

4、比较以下两种求+10x-2=0的根到三位小数所需的计算量： (1)在区间(0，1)内用二分法： (2)用迭代法=(2-)/10,取初值=0。

7用下列方法求f(x)=-3x-1=0在=2附近的根。根的准确性=1.87938524... (1)用Newton法； (2)用弦截法，取=2,=1.9; (3)用抛物线法，取=1,=3,=2。

解：(1)Newton法程序： clear

m=input('m='); x(1)=2;

f=inline('x^3-3*x-1'); f1=inline('3*x^2 - 3'); for n=1:m

x(n+1)=x(n)-f(x(n))/f1(x(n)); end x(n+1) 运行结果： m=10 ans =

1.879385241571817 (2)弦截法程序： clear

m=input('m='); x(1)=2; x(2)=1.9;

f=inline('x^3-3*x-1'); for n=2:m

x(n+1)=x(n)-(f(x(n))/(f(x(n))-f(x(n-1))))*(x(n)-x(n-1)); end x(n+1) 运行结果： m=9 ans =

1.879385241571817 (3)抛物线法程序： 首先建立M文件： function w=W(a,b) f=inline('x^3-3*x-1'); w=(f(b)-f(a))/(b-a); 再有程序： clear

m=input('m='); x(1)=1; x(2)=3; x(3)=2;

f=inline('x^3-3*x-1'); for n=3:m

q(n)=W(x(n),x(n-1));

p(n)=(W(x(n),x(n-2))-W(x(n),x(n-1)))/(x(n-2)-x(n-1)); w(n)=q(n)+p(n)*(x(n)-x(n-1));

s(n)=sqrt(w(n)^2-4*f(x(n))*p(n)); if w(n)+s(n)>0

x(n+1)=x(n)-2*f(x(n))/(w(n)+sqrt(w(n)^2-4*f(x(n))*p(n))); end

if w(n)-s(n)>0

x(n+1)=x(n)-2*f(x(n))/(w(n)-sqrt(w(n)^2-4*f(x(n))*p(n))); end end x=x(n+1) 运行结果： m=10 x =

1.879385241570596

13．应用Newton法于方程f(x)=1-=0，倒出求的迭代公式。 clear

m=input('m='); x(1)=11;

f=inline('1-115/x^2'); f1=inline('2*115*x^(-3)'); for n=1:m

x(n+1)=x(n)-f(x(n))/f1(x(n)); end

sqrta=x(n+1) 运行结果： m=10 sqrta =

10.723805294763608

实验报告七

一．实验名称

解线性方程组的直接方法

二．实验目的

掌握并运用直接法和迭代法

三．实验原理

线性方程组在应用数学中几乎处处都要遇到,因此线性方程组的数值解法必然成为计算数学最...为提高计算的精确度，适应计算机的特点，人们对古典消元法进行了一些重要的改进，这就形成了当前通用的方法——主元消去法，也叫高斯(Gauss)消元法

四． 实验内容

?4?1??1??1?0?已知矩阵

??An?????????14111114111114101114???????????????00000????????????????????????????????0?0000??0000??0000?0000????????????????????????000??00000???1?00000???1??00000???1?00000???0用最高精度的高斯分解，LU分解，平方根分解。

解： >> n=15 n = 15 >> a=ones(n);

b=diag(diag(a,-1),-1); c=diag(diag(a,-2),-2); d=diag(diag(a,1),1); e=diag(diag(a,2),2); f=diag(diag(a,-3),-3); g=diag(diag(a,3),3); h=4*eye(n); A=b+c+d+e+f+g+h A =

4 1 1 1 0 0 0 0 0 0 1 4 1 1 1 0 0 0 0 0 1 1 4 1 1 1 0 0 0 0 1 1 1 4 1 1 1 0 0 0 0 1 1 1 4 1 1 1 0 0 0 0 1 1 1 4 1 1 1 0 0 0 0 1 1 1 4 1 1 1 0 0 0 0 1 1 1 4 1 1 0 0 0 0 0 1 1 1 4 1 0 0 0 0 0 0 1 1 1 4

>> (1)LU分解： >> [L,U]=lu(A) L =

1.0000 0 0 0 0 0 0 0 0 0.2500 1.0000 0 0 0 0 0 0 0 411141114111 0 0

1?1??1?4??

0.2500 0.2000 1.0000 0 0 0 0 0 0 0 0.2500 0.2000 0.1667 1.0000 0 0 0 0 0 0 0 0.2667 0.2222 0.1905 1.0000 0 0 0 0 0 0 0 0.2778 0.2381 0.1806 1.0000 0 0 0 0 0 0 0 0.2857 0.2361 0.1805 1.0000 0 0 0 0 0 0 0 0.2917 0.2402 0.1805 1.0000 0 0 0 0 0 0 0 0.2931 0.2402 0.1799 1.0000 0 0 0 0 0 0 0 0.2931 0.2410 0.1801 1.0000 U =

4.0000 1.0000 1.0000 1.0000 0 0 0 0 0 0 0 3.7500 0.7500 0.7500 1.0000 0 0 0 0 0 0 0 3.6000 0.6000 0.8000 1.0000 0 0 0 0 0 0 0 3.5000 0.6667 0.8333 1.0000 0 0 0 0 0 0 0 3.4286 0.6190 0.8095 1.0000 0 0 0 0 0 0 0 3.4120 0.6157 0.8194 1.0000 0 0 0 0 0 0 0 3.4120 0.6160 0.8195 1.0000 0 0 0 0 0 0 0 3.4003 0.6119 0.8195 0 0 0 0 0 0 0 0 3.4000 0.6123 0 0 0 0 0 0 0 0 0 3.3991

(2)Cholesky分解： >> [R,p]=chol(A) R =

2.0000 0.5000 0.5000 0.5000 0 0 0 0 0 0 0 1.9365 0.3873 0.3873 0.5164 0 0 0 0 0 0 0 1.8974 0.3162 0.4216 0.5270 0 0 0 0 0 0 0 1.8708 0.3563 0.4454 0.5345 0 0 0 0 0 0 0 1.8516 0.3343 0.4372 0.5401 0 0 0 0 0 0 0 1.8472 0.3333 0.4436 0.5414 0 0 0 0 0 0 0 1.8472 0.3335 0.4437 0.5414 0 0 0 0 0 0 0 1.8440 0.3318 0.4444 0 0 0 0 0 0 0 0 1.8439 0.3321 0 0 0 0 0 0 0 0 0 1.8437 p = 0 (3)QR分解 >> [Q,R]=qr(A) Q =

-0.9177 0.2879 0.1837 0.1248 -0.1400 -0.0532 0.0072 0.0542 0.0080 0.0145 -0.2294 -0.9049 0.2391 0.1624 0.1488 -0.1297 -0.0546 -0.0023 0.0544 -0.0120 -0.2294 -0.1234 -0.8951 0.2240 0.1762 0.1576 -0.1309 -0.0675 -0.0075 -0.0559 -0.2294 -0.1234 -0.0787 -0.8855 0.2349 0.1851 0.1565 -0.1470 -0.0790 0.0100 0 -0.2605 -0.1662 -0.1129 -0.8663 0.2292 0.1854 0.1696 -0.1391 0.0796 0 0 -0.2721 -0.1848 -0.0822 -0.8618 0.2288 0.1956 0.1857 0.1315

0 0 0 -0.2773 -0.1761 -0.0824 -0.8621 0.2384 0.2146 -0.1976 0 0 0 0 -0.2896 -0.1858 -0.0794 -0.8956 0.1883 -0.1943 0 0 0 0 0 -0.2912 -0.1849 -0.0802 -0.9200 -0.1678 0 0 0 0 0 0 -0.2913 -0.1961 -0.0942 0.9316 R =

-4.3589 -2.2942 -2.2942 -2.2942 -0.6882 -0.4588 -0.2294 0 0 0 0 -3.8389 -1.4944 -1.4944 -2.1936 -0.5073 -0.3839 -0.2605 0 0 0 0 -3.6747 -1.2256 -1.6717 -2.2285 -0.5170 -0.4383 -0.2721 0 0 0 0 -3.6058 -1.4128 -1.7910 -2.2925 -0.5751 -0.4622 -0.2773 0 0 0 0 -3.4531 -1.2498 -1.7077 -2.2830 -0.5480 -0.4657 0 0 0 0 0 -3.4345 -1.2542 -1.7492 -2.2946 -0.5594 0 0 0 0 0 0 -3.4333 -1.2415 -1.7435 -2.2915 0 0 0 0 0 0 0 -3.2552 -0.9786 -1.5219 0 0 0 0 0 0 0 0 -3.1858 -0.8941 0 0 0 0 0 0 0 0 0 3.1665

实验报告八

一．实验名称

解线性方程组的迭代法

二． 实验目的

掌握并运用迭代法

三． 实验原理

雅可比迭代法是求解大型线性方程组的基本方法。利用GPU(Graphics Processing Unit,图形处理器)的并行处理能力,将雅可比迭代求解线性方程组过程中运算量较大的部分移植到GPU上执行,以提高运算速度。雅可比迭代法原理：

设方程组Ax=b 满足aii ≠0, 将方程组变形为: x=Bx+f, 则雅可比(Jacobi)迭代法是指 x

(k+1)

=Bx+f

(k)

即由初始解逐步迭代即可得到方程组的解。

四．实验内容

用Gauss-Seidel迭代法和Jacobi迭代法解线性方程组

Anx=b, 已知矩阵

?4?1??1??1?0???An?????????14111114111114101114???????????????00000????????????????????????????????0?0000??0000??0000?0000????????????????????????000??00000???14?00000???11??00000???11?00000???01 解：

(1)Jacobi迭代法

function [y,n]=jacobi(A,b,x0,eps) if nargin==10 eps=1.0e-6; elseif nargin<10 error return end

D=diag(diag(A)); L=-tril(A,-1); U=-triu(A,1); B=D\\(L+U); f=D\\b; y=B*x0+f; n=1;

while norm(y-x0)>=eps x0=y; y=B*x0+f; n=n+1; end

>> b=[0 0 0 0 0 0 0 0 0 0]'；

>> [x,n]=Jacobi(A,b,[0 0 0 0 0 0 0 0 0 0]',1.0e-6) x = 0 0 0 0 0 0

111?411??141?114??

0 0 0 0 n = 100

(2)Gauss-Seidel迭代法

function [y,n]=gauseidel(A,b,x0,eps) if nargin==10 eps=1.0e-6; elseif nargin<10 error return end

D=diag(diag(A)); L=-tril(A,-1); U=-triu(A,1); G=(D-L)\\U; f=(D-L)\\b; y=G*x0+f; n=100;

while norm(y-x0)>=eps x0=y; y=G*x0+f; n=n+1; end

>> b=[0 0 0 0 0 0 0 0 0 0]’;

>> [x,n]=jauseidel(A,b,[0 0 0 0 0 0 0 0 0 0]',1.0e-6) x = 0 0 0 0 0 0 0 0 0 0 n = 100

实验报告九

一． 实验名称

矩阵的特征植与特征向量计算。

二． 实验目的

掌握矩阵的特征植与特征向量计算方法。

三． 实验原理

线性代数中的矩阵运算；本实验所用Matlab命令提示：

(1)矩阵的输入格式：(2)求A的转置：A ’； (3)求A加B：A+B； (4)求A减B：A－B； (5)求数k乘以A：k*A； (6)求A乘以B：A*B； (7)求A的行列式：det(A)； (8)求A的秩：rank(A)； (9)求A的逆：inv(A)； (10)B右乘A的逆(11)B左乘A的逆

： B / A； ：A\\B；

；b=初始值 ：步长 ：终值；

(12)求A的特征值：eig(A)；

(13)求A的特征向量矩阵X及对角矩阵D：[X，D]=eig(A)； (14)方阵A的n次幂： A^ n ； (15)A与B的对应元素相乘：A .* B；

(16)存储工作空间变量：save’文件名’’变量名’； (17)列出工作空间所有变量：whos

四．实验内容

?4?1??1?已知矩阵

?1?0???An??????????0??0?0??014111????0000114111114101114???????????????0000000000?00000??00000???11410?0??0??0?0????????????1??1?1??4?

?????????0000?????????0000????????0000????????????111041111411????????求其最大特征值。

解：

[V,D]=eig(A,'nobalance') V =

0.3186 -0.0000 0.2581 0.3451 -0.4628 0.0182 -0.2318 0.5297 -0.3601 0.1852 -0.0478 -0.4082 -0.1944 -0.5734 -0.0675 -0.1813 0.3829 0.2620 -0.3859 0.2507 -0.4142 0.0000 0.2581 0.2283 0.5303 -0.3809 -0.2318 -0.0057 -0.3601 0.3162 -0.2211 0.4082 -0.4910 -0.0000 0.0000 0.5547 -0.1023 -0.0184 -0.2861 0.3782 0.4193 0.4082 0.2966 0.0000 0.0000 -0.1183 0.4852 -0.3879 -0.0998 0.3997 0.4193 -0.4082 -0.2966 -0.0000 -0.0000 -0.1183 -0.4852 -0.3879 0.0998 0.3997 -0.2211 -0.4082 0.4910 0.0000 0.0000 0.5547 0.1023 -0.0184 0.2861 0.3782 -0.4142 0.0000 -0.2581 0.3451 -0.4628 -0.3809 0.2318 -0.0057 0.3601 0.3162 -0.0478 0.4082 0.1944 -0.5734 -0.0675 -0.1813 -0.3829 0.2620 0.3859 0.2507 0.3186 -0.0000 -0.2581 0.2283 0.5303 0.0182 0.2318 0.5297 0.3601 0.1852 D =

1.8556 0 0 0 0 0 0 0 0 0 0 2.0000 0 0 0 0 0 0 0 0 0 0 2.3446 0 0 0 0 0 0 0 0 0 0 3.0000 0 0 0 0 0 0 0 0 0 0 3.0000 0 0 0 0 0 0 0 0 0 0 3.5929 0 0 0 0 0 0 0 0 0 0 3.7892 0 0 0 0 0 0 0 0 0 0 4.4490 0 0 0 0 0 0 0 0 0 0 6.8662 0 0 0 0 0 0 0 0 0 0 9.1024