数值分析_数学_自然科学_专业资料

1、一、算法设计方案1、解非线性方程组首先将x, y当作已知的常数，求解四个未知数t, u, w, uo利用newton法（简单迭 代法不收敛）求解非线性方程组，得到与x, y对应的向量t, uo求解步骤：1 ）> 选取初始向量t,u,v,w=l, 1 , 1 , 1；2）、计算咐>）和f'0）；3）、解关于心伙）的线性方程组（调用doolittle分解法求解此线性方程组）；4）、若 hil /|卅）ii",则取 x*"®；否则转 5；5）、计算兀"+d=jv（“）+ax；题冃中newton法迭代公式为：一 0.3sinx® 1

2、11时厂0.5cos兀$）+兀2“）+兀3+兀伙）x- 2.6710.5 cos 勺伙)11时）护）+0.5sinx2 +屮）+兀4®-丹_1.070.51 -sin 寸)1时）0.5西 +兀2伙)4-cosx3a,+ x4(a) -x.-3.7410.51cosx4a)xyk) +0.5x2(a>+x3m) +sinx4(a)-yj-0.79口屮坷=0.08z, yj = 0.5 + 0.05 j（i = 0,1,2,10;丿=0,1,2，. 20）2、分片二次代数插值解题思路：由1得到的x, y和t, u的映射表，f（t（x, y）, u（x, y）,即求得f （x, y）

3、。但由于得到 的t, u不可能正好是题日提供的二维数表屮的值，需耍用相关规则对插值节点加以规范。利用（x, y）以及对应的f（x, y）,就可能通过二元拉格朗li插值多项式得到f（x, y） 的表达式。插值节点:1）、根据计算得到的t、u值，选取插值节点；选择标准如下：假设对（t, u）,这里用（x, y）代替：设：x： =4- ih i = 0,1,2,.,丹='o + 丿=°丄2,，加a）、若满足：xi-h/2<x<xj-h/2<i<n-yj-t/2<y<yj-t/2i< j<m-则应选择（无，片）伙='_ 口丿+1

4、,厂二丿° _ 1, j, j +1）为插值节点b）、若满足:x<x -0.5hx>+°.5，则取= 1 或,=/7_1:y < 刃一0.51或 y > x-i + 0.5丁，则取丿=1 或丿=加_ 1 ；2）、双元二次插值子程序相应的插值多项式为：022（九刃二 x £人（必（刃/（母,儿）k=il r=j-l其中山+1tk3、最小二乘法曲血拟合设在三维直角坐标系oxyu中给定（ni+l） *（n+l）个点（即三维坐标）(4,儿,s)(z = 0,l,.,m;j = 0,l,.,h)在本题中即为（"力"（*"

5、）。选定m+1个x的函数0（兀）（厂=0丄,"）以及n+1 个y的函数0$o）g = °丄,"）。本题中卩（力二*,肖$（刃=讯，n = m=kf于是得到乘积型基函数构成的曲面，kp（兀,刃=为c/yr,.v=010 20g m /a, y）- pg, x）2随着k值的不断增大，粘度匸0戶0会越来越大，题冃要求精度为10_7,此时的k即为要求的最小值。解题思路：1）、求解矩阵a固定力，以©（x）二"为基函数对数据（*知）作最小二乘拟合，得到n+1条拟合曲线r=0其中（购八夠,r=勺是法方程bt baj = btuj j = 0,1,.,h的解，

6、而"二©3）（“zx（如），求解n+1线性方程组，得到矩阵a。2）、求解矩阵gg 0 （刀）（“+）x（n+l）3）、系数矩阵cc = a（gtgy'gtr4、子程序说明了程序名称功能subroutine f fit(tl, t2, c, sigma)最小二乘法曲而拟合子程序，可给出拟合精度sigmasubroutine f pxy (c, tl, t2, x, y, p xy)以x,y的幕函数为基，得到拟合系数矩阵csubroutine f zxy(z)分片插值子函数,利用已知的（x,y），得到z（x,y）subroutine f zut(u, t, p)分片插值

7、子函数,利用求取的（u,t）,得到z（u,t）subroutine dlu (a, b, x)doolittle分解求线性方程组子函数subroutine f newton iteration(x,y, u,t)newton迭代法解非线性方程组子程序5、主程序main功能说明主程序对xi,yi赋值,通过调用子程序对非线性方程组求解,得到相应的数据(t,u,v,w), 通过调用插值了程序，得到对丿应的z=f(x,y),并以文件的形式进行输出。通过调用拟合了程序 对拟合系数短阵及拟合精度的求解，结果以文件形式输出。二、fortran源程序! / illi面拟合子两数，并给出拟合精度/subrout

8、ine f_fit(tl, t2, c, sigma)use imslimplicit noneinteger i, j, 11, t2parameter nl=llparameter n2=21dimension b(nl, 11), b_trans(t 1, nl), b_trans_b(tl, tl), b_inverse(tl, 11)dimension g(n2, t2), g trans (t2, n2), g trans g(t2, t2), g inverse(12, t2)double precision b, g, b trans,g trans, b trans_b, g

9、_trans_g, b_inverse, g_inverse dimension temp i(tl, nl), temp_2(tl, n2), temp_3(tl, t2), &c (tl, t2), u(nl, n2), p_xy (nl, n2), x(nl), y (n2)double precision temp_l, temp_2, temp_3, c, u, p_xy, sigma, x, y!/初始化x, y/do i=l,nlx(i)=0. 08*(i-l)end dodo j=l, n2y(j)=o. 5+0. 05*(j-1)end do!/据题意，求出矩阵b,

10、g/do i=l, nldo j 二 1,11b(i, j)=x(i)*(j-l)end doend dodo i=l, n2do j=l, t2g(i, j)=y(i)*(j-l)end doend do!/确定b, g的转置b_trans和g_trans/do i=l, nldo j=l,tlb_trans(j, i)=b(i, j)end doend dodo i=l,n2do j=l, t2g_trans(j, i)=g(i, j)end doend do!/求解 b_trans_b 和 g_treins_g 的逆矩阵 binverse 和 g_inverse/ b_trans_b=m

11、atmu 1 (b_trans, b) imatmul 为矩阵相乘函数，库函数 g_trans_g=matmul(g trans, g)b_inverse=. i. b_trans_bg_inverse二.i. g_trans_g!/求解 c 矩阵/call f_zxy(u)temp_l=matmul(b_inverse, b trans)temp_2=matmul(temp_l, u) temp_3=matmul(temp_2,g) c=matmui(temp 3,g inverse)!/求拟合精度误差/sigma=0do i=l, nldo j=l,n2call f_pxy(c, tl,

12、t2, x(i), y(j), p_xy(i, j)sigma=sigma+(u(i, j)p_xy(i, j)*2end doend dowrite (*, *) sigmaend subroutine f_fit!/子函数f_pxy给定拟合曲面的近似表达式/subroutine f_pxy (c, 11, t2, x, y, p_xy)implicit noneinteger 11, t2dimension c(tl, t2)double precision p_xy, temp_4, x, y, cinteger i, j, ktemp 4=0. ododo i=l,tldo j=l,

13、t2temp_4=temp_4+c (i, j) *(x*(it)*(y* (j-1)end doend dop_xy=temp_4 !拟合曲面end subroutine f_pxy!/分片插值子函数,得thz=z(x, y) /subroutine f_zxy(z)dimension x(11), y(21), u(ll, 21), t (11, 21), z(ll, 21) double precision x, y, u, t, zinteger i, jdo i = l, 11do j 二 1,21x(i)二0. 08*(i-l)y(j)=o. 5+0. 05* (j-1)!调用牛顿

14、迭代法求非线性方程子函数call f_newton_iteration(x(i), y(j), u(i, j), t (i, j) call f_zut(u(i, j), t(i, j), z(i, j)end doend doend subroutine f_zxy!/ doolittle分解子函数求解线性方程组/ subroutine dlu(a, b, x)integer, parameter : kkk=4real (kind二8) mreal (kind=8), dimension(kkk, kkk), intent(in) : areal(kind二8), dimension(kk

15、k, kkk):1=0,u=0real(kind=8), dimension(kkk): yreal(ki nd=8), dimension (kkk), intent (out):xreal (kind二8), dimension(kkk), intent(in):binteger i, j, kdo k=l, kkk, 1do j=k, kkk, 1m=0do t=l, k-1, 1m=l (k, t)*u(t, j)+mend dou(k, j)=a(k, j)-mend dodo i 二k+1, kkk, 1m=0do t=l, kt, 1m=m+l (i, t)*u(t, k)end

16、 do1 (i, k) = (a(i, k) -m) /u (k, k)end dol(k,k)=lend do!/解方程/y (l)=b(l)do i=l, kkk, 1m=0do t=l, i-1, 1m=m+l(i, t)*y(t)end doy (i)=b(i)-mend dox(kkk)=y(kkk)/u (kkk, kkk)do i=kkk-l, 1, -1m=0do t=i + l, kkk, 1in二m+u(i, t)*x (t)end dox(i) = (y(i)-m)/u(i, i)end doend subroutine dlu!/牛顿迭代法子程序/ subroutine

17、 f_newton_iteration(x, y, u, t) parameter n二4dimension f (n), aa(n, n), f delta(n) double precision f, aa, f delta double precision t, u, v, w, x, y double precision epsion, sl,s2integer i, j!迭代初始值t=l. 0u=1.0v=l. 0w=1.0epsion二1.0do while (epsion. ge. le12)f (l)=-(0. 5*cos(t)+u+v+w-x-2. 67)f (2)=-(t+

18、0. 5*sin (u) +v+w-y-l. 07)f (3)=-(0. 5*t+u+cos (v) +w-x3. 74)f(4)=-(t+0. 5*u+v+sin(w)-y-0. 79)aa(l, l)=-0. 5*sin(t)aa(l, 2)=1aa(l, 3) = 1aa(l,4)=laa(2, 1)=1aa(2, 2)=0. 5*cos (u)aa(2, 3) = 1aa 4) = 1aa(3, 1)=0. 5aa(3, 2)=1aa(3, 3)=-sin (v)aa(3, 4)=1aa(4, 1) = 1aa(4, 2) =0. 5aa(4, 3)=1aa(4, 4) =cos (w

19、)call dlu(aa, f, f_delta)sl=f_delta(l)*f_delta(l)+f_delta(2)*f_delta(2) +f_delta(3)*f_delta(3) +f_delta(4)*f_delta(4)sl=sqrt (si)s2二t*t+u*u+v*v+ws2=sqrt(s2)epsion二sl/s2t=t+f delta(l)u=u+f delta(2)v=v+f_delta(3)w=w+f_delta(4)end doend subroutine f newton iteration!/分片插值子函数,得出z=z(u, t)/ subroutine f_z

20、ut(u, t, p)implicit noneparameter n=6dimension x(n),y(n)double precision x,y, h, hldouble precision u, tinteger i, j, ii, jj!ii, jj 插值节点integer k, r, tidimension lx(n), ly(n) ,z(n, n)double precision p, lx, ly, z !p, lx, ly, z 插值多项式屮的项 x(l)=0y (1)=0h=0. 4hl=0. 2do i二2, nx(i)=x (l) + (il)*hy(i)=y(l)

21、+ (i-l)*hlend doii=0jj=op二 0!/给定z(n, n)/z(:, l) = (/-0. 5,-0. 34, 0. 14, 0. 94, 2. 06, 3. 5/)z(:,2) = (/-0. 42, -0. 5,-0. 26, 0. 3, 1. 18, 2. 38/)z(:,3) = (/-0. 18, -0. 5, -0. 5,-0. 18, 0. 46, 1.42/)z(:,4) = (/0. 22, -0. 34,-0. 58, -0. 5, -0. 1, 0. 62/)z(:,5) = (/0. 78,-0. 02, -0. 5, -0. 66, -0. 5,-

22、0. 02/)z(:,6) = (/1.5, 0. 46, -0. 26,-0. 66,-0. 74, -0. 5/)do i=3, n-2if (u. gt. (x(i)-h/2). and. (u. le. (x(i)+h/2) then 11=1goto 10else if (u. le. (x(3)-h/2) thenii=2elseii=n-lend ifend do10 do j=3, n-2if (t. gt. (y(j)_hl/2). and. (t. le. (y(j)+hl/2) then jj=jgoto 20else if (t. le. (y(3)-hl/2) the

23、njj=2elsejj=n-lend辻end do20 do k二ii-1, ii+1do r=jj-l, jj+1lx(k)=lly(r)=ldo ti=ii-l, ii+1if (ti. ne. k) thenlx (k) =lx (k) * (ux (ti) /(x (k) -x (ti)end ifend dodo ti=jj-l, jj+1if (ti. ne. r) thenly (r) =ly (r) * (t-y (ti)/(y (r) -y (ti)end ifend dop二p+lx (k)*ly (r) *z (k, r)end doend doend subroutin

24、e fzut!/主函数 main/program mainimplicit nonodouble precision vector x(0:10), vector_y(0:20)double precision xx(0:7),yy(0:4)double precisiontu (0:230, 0:l),ut(0:l),zt(0:230), c(0:10, 0:10), tu2 (45, 0:1), zt2(0:39), tt (45) double precision z,sigmainteger i, j, in, n, 1, t2, n2n=10m 二 20do i=0, 10vecto

25、r_x(i)=0. 08*iend dodo j=0, 20vector_y(j)=0. 5+0. 05*jend doopen (11, f 订e二'xytuz. txt')write (11, 100)' x',' y',' t',' u',' z'do i二0, 10do j=0, 20call f_newton_iteration(vector_x(i), vector_y(j), &close(33)t2h8f12h5doh-lrt-2xx(illvp1*end dodu n2yy

26、(l)h0 5+p2沃 jend dodoh-lrt-2do jhl n2call fnewtoniberation(xx(il). yy(jl)" tu2(il)*5+j"ptu2(i 丄¥5+jo)end doend dodoh-lkocall fzut(tu2(i-ptu2p0)9 zb2(ill)end dodo71-1t2do j*lyn2call fpxyf6" 6xx(i'l yy(j丄)"tt (l-l)*5+j)tu (i 決21+jyl)，tu (i 決21+jo)call fzut(tu(i*21+jy1)9tu(

27、i*21+j- 0)- zt (i*21+j) wrice(il 200) vecborx(i)" vectory(j)"gtu (黨21+jyptu (i 关21+y 0l zt (721+j)end doend do close(ll)open (22. f i le'lf (xy) bxt “)write (22400)" xi"yi"f (xr-yi)-dotrpndo jupmwrice (22- 300) vectorlx (i)，vectory (jl zt (*21+j)end doend doclose (22)do

28、1h16call f if i tpjpsi gma)end doopen (33- f ilc'lmatrixic, txt")doh-p5write (33.(20e20 12)9) (c (i, jl j6 5)end doend doopen(44, file二'f2(x, y). txt，)write(44, 600)' x*i',' y*i',' f(xi, yi)',' p(xi, yi)'do i=l, t2do j=l, n2write (44, 500) xx (i-1), yy (j

29、l), zt2 (it) *5+jt), tt (il) *5+j) end doend do close (44)! /输出格式控制/100 format(2al0, 3al5)200 format (2el2. 4, 3e20. 12)300 format (2el2.4, e20. 12)400 format(2al0, al8)500 format (2e20. 7, 2e20. 12) 600 format(2al8, 2al8) end program main 三、结果汇总1、数表:xi, yi、f (xi, yi)xi0.0000e+000.0000e+000.0000e+00

30、0.0000e+000. 0000e+00yi0.5000e+000.5500e+000.6000e+000.6500e+000.7000e+00f(xi,yi)0. 446504069241e+000.324683310597e+000.210159730631e+000.103043643055e+000.340192524000e-020.0000e+000.7500e+00 -0.887357886680e-010.0000e+000. 8000e+00 -0.173371611924e+000. 0000e+000.8500e+00 -0 250534594048e+000. 000

31、0e+000.9000e+00 -0. 320276491493e+000. 0000e+000.9500e+00 -0. 382668052743e+000. 0000e+000.1000e+01-0.437795745638e+000.0000e+000. 1050e+01-0.485758921243e+000.0000ei000. 1100ew10. 526667236102ew00. 0000e+000.1150e+01-0.560638463024e+000.0000e+000. 1200e+01-0.587796524654e+000.0000e+000. 1250e+01&qu

32、ot;0.608269768328e+000.0000e+000.1300e+01-0.622189446194e+000.0000e+000. 1350e+01-0.629688384281e+000.0000e+000. 1400e+01-0.630899769389e+000.0000e+000. 1450e+01-0.625956164380e+000.0000e+000. 1500e+01"0.614988550047e+000.8000e-010.5000e+00 0.638015234587e+00 0.8000e-01 0.5500e+00 0.50661179587

33、6e+000. 8000e-01 0. 6000e+00 0. 382176408075e+008000e-018000e018000e-018000e-018000e-018000e-018000e018000e-018000e-018000e-018000e-018000e-018000e-018000e-018000e-018000e-018000e-018000e-010.6500e+00 0. 264863525560e+000. 7000e4 00 0. 154780230633e-1000. 7500e+00 0. 519927097057e-010. 8000e+00 -0.

34、434680179481e-010.8500e+00 -0. 131601038093e+000.9000e+00 -o. 212431072587e+000- 9500e+00 -0.286004537561e+000.1000e+010.1050e+010.1100e+010. 1150e+010.1200e+010. 1250e+010.1300e+010. 1350e+010.1400e+01-0.352386059629e+00-0.411655437858e+00"0.463904893928e+00-0.509236708669e+00-0.547761104059e+

35、00-0.579594377202e+00-0.604857251213e+00-0.623673426166e+00-0.636168257103e+000. 1450e+01 -0.642467668455e+000-1500e+01 -0. 642697116743e+001600e+000. 5000e+001600e+000. 5500e+001600e+000. 6000e+001600e+000.6500e+001600e+000.7000e+001600e+000.7500e+001600e+000.8000e+001600e+000. 8500e+000. 840081396

36、883e+000.699764166271e+000. 566061444836e+000. 439171614444e+000. 319242167288e+000.206376218657e+000.100638546919e+000.206075984654e-021600e+00 0.9000e+00 -0. 893540080136e-011600e+001600e+001600e+001600e+001600e+001600ei001600e+001600e+001600e+001600e+001600e+001600e+002400e+002400e+002400e+002400

37、e+000.9500e+00 -0. 173626954709e+000.1000e+010.1050e+010.1100e+010. 1150e+010-1200e+010.1250e+010. 1300e+010. 1350e+010.1400e+010. 1450e+010. 1500e+010.5000e+000.5500e+000.6000e+000.6500e+00-0.250799943998e+00-0.320932252750e+00-0.384097719188e+00-0.440382160815e+000.489881139368e»00-0.53269795

38、4784e+00-0.568941871351e+00"0.598726545062e+00-0.622168637157e+00-0.639386546614e+00-0.650499364497e+000.105151509248e+010.902927427685e+000.760580262813e+000.624715195878e+002400e+002400e4002400e+002400e+002400e+002400e+002400e4002400e+002400e+002400e+002400e+002400e+002400e+002400e+002400e+00

39、2400e+003200e+003200e+003200e+003200e+003200e+003200e+003200e+003200e+003200e+003200e+003200e+003200e+003200e+003200e+003200e+003200e+003200e4003200e+003200e+003200e+003200e+004000e+004000e+004000e+004000e+004000e+004000e+004000e+000.7500e+000. 8000ew00.8500e+000.9000e+000.373134067685e+000. 2576567

40、76698e+000.149150588766e+000.476470043700e-010.9500e+00 -0.468493081747e-010.1000e+010-1050e+010. 1100e+010. 1150e+010.1200e+010. 1250e+010. 1300e+010. 1350e+010.1400e+010.1450e+010.1500e+010.5000e+000- 5500e+000.6000e+000.6500e+000.7000e+000.7500e+000.8000e+000-8500e+000.9000e+000.9500e+000.1000e+0

41、10.1050e+010. 1100e+010. 1150e+010.1200e+010. 1250e+010-1300e+010. 1350e+010. 1400e+010. 1450e+010.1500e+010.5000e+000- 5500e+000.6000e+000.6500e+000.7000e+000. 7500e+000-8000e+00-0.134356747670e+000.214913334079ei00-0.288573687009e+00-0.355406352189e+00"0.415491385199e+00-0.468918240380e+00-0.

42、515783875727e+00-0.556191070573e+00-0.590246929357e+00-0.618061557834e+00-0.639746851961e+000.127124675843e+010. 111500201761e+010.964607722308e+000.820347363242e+000.682447673529e+000.551085227171e+000.426392408546e+000.308463021403e+000.197357157387e+000.931056491500e-01-0.428596469251e-02-0.94833

43、9139383e-01-0.178572979137e+00-0.255553768875e+00-0.325840141101e+00-0.389506980643e+000.446638198547e+00-0.497324947822e+00-0.541664031390e+00"0.579756487909e+00-0.611706299458e+000.149832107231e+010.133499864071e+010. 117712512422e+010.102502405047e+010.878960016740e+000.739145103470e+000.605

44、744887693e+004000e+000.8500e+004000e4000. 9000e i004000e+000.9500e+004000e+000.1000e+014000e+000. 1050e+014000e+000. 1100e+014000e4000. 1150e+014000e+000. 1200e+014000e+000. 1250e+014000e+000. 1300e+014000e+000. 1350e+014000e+000.1400e+014000e+000.1450e+014000e+000.1500e+010. 478883880984e+000. 3586

45、50648829e+000.245102261159e+000.138268376673e+000.381548905914e-01 -0.552527979027e-010.141986870438e+00 -0.222094430862e+00 -0.295635227198e+00 -0.362679507525e+00 -0.423306161793e+00 -0.477601035622e+00 -0.525655428548e+00 -0.567564753485e+004800e+004800e+004800e+004800e+004800e+004800e+004800e+00

46、4800e+004800e+004800e+004800e+004800e+004800e+004800e+004800e+004800e+004800e+004800e+004800e4004800e+004800e+005600e+005600e+005600e+005600e+005600e+005600e+005600e+005600e+005600e+000. 5000e+000. 5500e+000. 6000e+000. 6500e+000. 7000e+000. 7500e+000. 8000e+000.8500e+000.9000e+000.9500e+000.1000e+0

47、10.1050e+010.1100e+010. 1150e+010. 1200e+010. 1250e+010. 1300e+010. 1350e+010. 1400e+010.1450e+010.1500e+010.5000e+000.5500e+000.6000e+000.6500e+000.7000e+000.7500e+000.8000e+000.8500e+000.9000e+000.173189277935e+010.156203460147e+010.139721693052e+010.123780101007e+010.108408753012e+010.93632276707

48、9e+000.794704459103e+000.659387211130e+000.530487603091e+000.408088704433e+000.292244222111e+000.182982228522e+000.803085151786e-01 -0.157904052160e-01 -0.105344546036e+00 -0.188398086861e+00 -0.265007147050e+00 -0.335237837830e+00 -0.399164503807e+00 -0.456868144112e+00 -0.508435001840e+000.1971221

49、85031e+010.179532964555e+010.162406714373e+010.145783060027e+010.129695465723e+010.114171810591e+010.992349529967e+000.849032670029e+000.711911361038e+005600e+005600e4005600e+005600e+005600e+005600e+005600e4005600e+005600e+005600e+005600e+005600e+006400e+006400e+006400e+006400e+006400e+006400e+00640

50、0e+006400e+006400e+006400e+006400e+006400e+006400e+006400e+006400e+006400e+006400e+006400e+006400e+006400e+006400e4007200e+007200e+007200e+007200e+007200e+007200e+007200e+007200e+007200e+007200e+007200e+000.9500e+000. 1000ew10.1050e+010.1100e+010. 1150e+010. 1200e+010- 1250e+010. 1300e+010. 1350e+01

51、0.1400e+010. 1450e+010.1500e+010-5000e+000.5500e+000.6000e+000.6500e+000.7000e+000-7500e+000. 8000e+000. 8500e+000.9000e+000.9500e+000.1000e+010.1050e+010. 1100e+010. 1150e+010.1200e+010. 1250e+010.1300e+010. 1350e+010.1400e+010. 1450e+010-1500e+010.5000e+000.5500e+000.6000e+000.6500e+000.7000e+000-7500e+000.8000e+000.8500e+000.9000e+000.9500e+000.1000e+010. 581094170177e+000. 456658526807e+000.338654511410e+000.227108271757e+000.122025104649e+000.233922218851e-01-0.688187015106e-01-0.154649345159e+00-0.234152668459e+00"0.307391094659e+00-0 374434865587e+00-0.4353