MATLAB 程序1、图示牛顿迭代法（M 文件）文件名：newt_gfunction x = new_g(f_name,x0,xmin,xmax,n_points)clf,hold off% newton_method with graphic illustrationdel_x = 0.001;wid_x = xmax - xmin; dx = (xmax - xmin)/n_points;xp = xmin:dx:xmax; yp = feval(f_name,xp);plot(xp,yp);xlabel('x');ylabel('f(x)');title('newton iteration'),hold onymin = min(yp); ymax = max(yp); wid_y = ymax-ymin;yp = 0. * xp; plot(xp,yp)x = x0; xb = x+999; n=0;while abs(x-xb) > 0.000001if n > 300 break; endy=feval(f_name,x); plot([x,x],[y,0]);plot(x,0,'o')fprintf(' n = % 3.0f, x = % 12.5e, y = % 12.5e \ n', n, x, y);xsc = (x-xmin)/wid_x;if n < 4, text(x,wid_y/20,[num2str(n)]), endy_driv = (feval(f_name,x + del_x) - y)/del_x;xb = x;x = xb - y/y_driv; n = n+1;plot([xb,x],[y,0])endplot([x x],[0.05 * wid_y 0.2 * wid_y])text( x, 0.2 * wid_y, 'final solution')plot([ x ( x - wid_x * 0.004)], [0.01 * wid_y 0.09 * wid_y])plot([ x ( x + wid_x * 0.004)], [0.01 * wid_y 0.09 * wid_y])传热问题假设一个火炉是用厚度为0.05m 的砖单层砌成的。

炉内壁温度为T 0=625K, 外壁温度为T 1（未知）。

由于对流和辐射造成了外壁的热量损失，温度T 1由下式决定：44111()()()()0f k f T T T T T h T T xεσ∞=-+-+-=∆ 其中：k ：炉壁的热传导系数，1.2W/mKε: 发射率，0.8T 0：内壁温度，625KT 1：外壁温度（未知），KT ∞：环境温度，298KT f ：空气温度，298KH ：热交换系数，20W/m 2Kσ:Stefan-Boltzmann 常数，5.67*10-8W/m 2K 4x ∆:炉壁厚度，0.05m试用牛顿迭代法求解T 1.先根据待求解的方程式建立（M 文件）文件名：wall_htfunction f= wall_ht(t1)k = 1.2; e = 0.8; tinf = 298;tf = 298; h = 20; t0 = 625;sig = 5.67e-8; wall_thick = 0.05;f = k/wall_thick * (t1-t0) + e * sig * (t1.^4-tinf^4) +h * (t1 - tf);程序的调用及计算结果：>> newt_g('wall_ht',550,400,600,500)n = 0, x = 5.50000e+002, y = 7.03301e+003n = 1, x = 4.55199e+002, y = 6.58551e+002n = 2, x = 4.44423e+002, y = 6.44623e+000n = 3, x = 4.44316e+002, y = 6.27680e-004n = 4, x = 4.44316e+002, y = 5.67525e-010ans =444.3157平衡问题一氧化碳与氢按以下反应生成甲醇232fr K K CO H CH OH + 现有1molCO 与2molH 2的混合物，在温度t=5900C ，压力p=3.04*107Pa 下进行反应并达到平衡（Kf=1.393*10-15，Kr=0.43）,求CH3OH 在平衡气中的摩尔分数。

达到平衡时的关系式：232222()4(1)12(1)(129)40f x Kp x Kp x Kp x Kp =+-+++-=其中fr K K K据待求函数建立（M 文件）文件名：ping_henfunction f= ping_heng(x)t = 590; p = 3.04e7; kf = 1.393e-15;kr = 0.43;f = 4 * ((kf / kr) * p^2 + 1) * x.^3 - 12 * ((kf / kr) * p^2 + 1) * x.^2 + (12 * (kf / kr) * p^2 + 9) * x - 4 * (kf / kr) * p^2程序的调用及计算结果：newt_g('ping_hen',0.5,0.0,1.0,500)f = 0.5031n = 0, x = 5.00000e-001, y = 5.03076e-001f = 0.5120f = -0.0798n = 1, x = 4.43838e-001, y = -7.97582e-002f = -0.0680f = -0.0010n = 2, x = 4.50599e-001, y = -1.03321e-003f = 0.0104f = 2.1619e-006n = 3, x = 4.50689e-001, y = 2.16189e-006f = 0.0114ans =0.45072、线性方程组求解>> a=[3 2; 1 -1]a =3 21 -1>> y = [-1, 1]';>> x=a\yx =0.2000-0.8000高斯消去法求解线性方程组-0.0400x1 + 0.0400x2 + 0.1200x3= 3.00000.5600x1 - 1.5600x2 + 0.3200 x3= 1.0000-0.2400x1 + 1.2400x2- 0.2800 x3 = 0>> clear>> a = [-0.04 0.04 0.12 3; 0.56 -1.56 0.32 1; -0.24 1.24 -0.28 0] a =-0.0400 0.0400 0.1200 3.00000.5600 -1.5600 0.3200 1.0000-0.2400 1.2400 -0.2800 0>> x = [0,0,0]'x =>> tempo = a(1,:);a(1,:) = a(2,:);a(2,:) = tempo;aa =0.5600 -1.5600 0.3200 1.0000-0.0400 0.0400 0.1200 3.0000-0.2400 1.2400 -0.2800 0>>>> a(2,:) = a(2,:) - a(1,:) * a(2,1)/a(1,1);>> a(3,:) = a(3,:) - a(1,:) * a(3,1)/a(1,1);aa =0.5600 -1.5600 0.3200 1.00000 -0.0714 0.1429 3.07140.0000 0.5714 -0.1429 0.4286>> tempo = a(3,:);a(3,:) = a(2,:);a(2,:) = tempo;aa =0.5600 -1.5600 0.3200 1.00000.0000 0.5714 -0.1429 0.42860 -0.0714 0.1429 3.0714>>a(3,:) = a(3,:) - a(2,:) * a(3,2)/a(2,2);aa =0.5600 -1.5600 0.3200 1.0000 0.0000 0.5714 -0.1429 0.4286 0.0000 0 0.1250 3.1250 >> x(3) = a(3,4)/a(3,3);>> x(2) = (a(2,4) - a(2,3) * x(3))/a(2,2);>> x(1) = (a(1,4) - a(1,2:3) * x(2:3))/a(1,1);xx =7.00007.000025.00003、数值积分梯形法（M 文件）文件名：trapez_v% function trapez_v(f,h) integrates a function defined in % vector f with interval size, h.% Copyright S. Nakamura, 1995function I = trapez_v(f, h)I = h*(sum(f) - (f(1) + f(length(f)))/2);对方程220008.11200du u u dx=--采用梯形法积分 >> clear>> n_points = 16;i = 1:n_points;>> h = (30-15)/(n_points-1); u = 15 + (i-1) * h; >> f = 2000 * u./(8.1 * u.^2 + 1200);>> i = trapez_v(f,h)i =127.5040。