上海交通大学计算方法大作业
计算方法上机报告_董晓壮
计算方法A上机报告学院（系）：电气工程学院学生姓名：陶然学号：授课老师：完成日期：2019年12月03日西安交通大学Xi'an Jiaotong University目录1 QR分解法求解线性方程组 (2)1.1 算法原理 (2)1.1.1 基于吉文斯变换的QR分解 (2)1.1.2 基于豪斯霍尔德变换的QR分解 (3)1.2 程序流程图 (4)1.2.1 基于吉文斯变换的QR分解流程图 (4)1.2.2 基于豪斯霍尔德变换的QR分解流程图 (5)1.3 程序使用说明 (5)1.3.1 基于吉文斯变换的QR分解程序说明 (5)1.3.2 基于豪斯霍尔德变换的QR分解程序说明 (7)1.4 算例计算结果 (8)2 共轭梯度法求解线性方程组 (10)2.1 算法原理 (10)2.2 程序流程图 (10)2.3 程序使用说明 (11)2.4 算例计算结果 (12)3 三次样条插值 (14)3.1 算法原理 (14)3.2 程序流程图 (16)3.3 程序使用说明 (17)3.4 算例计算结果 (19)4 龙贝格积分 (21)4.1 算法原理 (21)4.2 程序流程图 (22)4.3 程序使用说明 (23)4.4 算例计算结果 (24)结论 (26)1 QR 分解法求解线性方程组1.1 算法原理矩阵的QR 分解是指，可以将矩阵A 分解为一个正交矩阵Q 和一个上三角矩阵R 的乘积，实际中，QR 分解经常被用来解决线性最小二乘问题，分解情况如图1.1所示。
=⨯图1.1 QR 分解示意图本次上机学习主要进行了两个最基本的正交变换—吉文斯（Givens ）变换和豪斯霍尔德（Householder ）变换，并由此导出矩阵的QR 分解以及求解线性方程组的的方法。
1.1.1 基于吉文斯变换的QR 分解吉文斯矩阵也称初等旋转阵，如式（1.1）所示，它把n 阶单位矩阵I 的第,i j 行的对角元改为c ，将第i 行第j 列的元素改为s ，第j 行第i 列的元素改为s −后形成的矩阵。
计算方法第7章／《数值分析》／清华大学／上海交通大学／西安交通大学
Euler-Maclaurin 公式
å ò Tn - I =
k
c2 j [ f (2 j-1) (b) - f (2 j-1) (a)]h2 j +
b a
P2k ( x) f (2k ) (x)dx
j =1
8
å c2 j
=
(-1)
j +1
1
(2p )2
j
¥ k =1
1 k2j
，|
P2k (x) |£ c2k h2k
同上
ò RS =
b f (4) (x ) (x - a)(x - c)2 (x - b)dx a 4!
=- b - a (b - a )4 f (4) (h) 180 2
5
复化公式及误差分析
由上述误差表达式可知，区间越小，绝对误差越小，复化梯形公式：
将积分区间
n
等分，节点是
xi
=
a
+ ih, h
=
值公式
pn (xk ) = f (xk ) = yk
利用 Lagrange 插
1
å Õ pn (x) =
nn
(
k=0 j=0 j¹k
x xk
- xj - xj
)yk
¬ 做代换 x = a + th,t
=
x-a
å Õ n n t - j
=
(
k =0
j=0
k
-
j )yk
h
j¹k
以 pn (x) 代 f (x) 得
k -1
<e
停止
输出 Tk(k )
»
I
。否则 h
Ü
h 2
计算方法大作业
计算方法大作业学院：电子工程姓名：班级：学号：大作业选题：分析方程求根问题中牛顿法的性能，包括收敛性等，并用该方法求解一个问题，给出过程和结果。
一、牛顿迭代法介绍：用迭代法求方程0)(=x f 的根时，首先要构造一个迭代函数，迭代函数构造的好坏，不仅影响收敛速度，而且有可能使迭代序列发散，构造迭代函数的一条重要途径，是用近似方程代替原方程去求根，因此如果能将非线性方程0)(=x f 用线性方程来近似代替，那么求近似根问题就容易得到解决，而且十分方便。
牛顿法就是把非线性方程线性化的一种方法。
二、牛顿迭代法原理设已知方程0)(=x f 的近似根0x ,则在0x 附近)(x f 可用一阶泰勒多项式))((')()(000x x x f x f x p -+=近似代替.因此, 方程0)(=x f 可近似地表示为0)(=x p .用1x 表示0)(=x p 的根,它与0)(=x f 的根差异不大.设0)('0≠x f ,由于1x 满足,0))((')(0100=-+x x x f x f 解得)(')(0001x f x f x x -= 重复这一过程,得到迭代公式)(')(1n n n n x f x f x x -=+ 这就是著名的牛顿迭代公式,它相应的不动点方程为)(')()(x f x f x x g -=.用牛顿迭代公式求方程根的方法称为牛顿迭代法，简称牛顿法。
三、牛顿迭代法的几何解析在0x 处作曲线的切线，切线方程为))((')(000x x x f x f y -+=。
令0=y ，可得切线与x 轴的交点坐标)(')(0001x f x f x x -=，这就是牛顿法的迭代公式。
因此，牛顿法又称“切线法”，其几何意义即为0x 点处的切线方程。
四、牛顿迭代法的收敛性 计算可得2)]('[)(")()('x f x f x f x g -=,设*x 是0)(=x f 的单根,有0)(',0)(**≠=x f x f ,则0)]('[)(")()('2****=-=x f x f x f x g , 故在*x 附近,有1)('<x g .根据不动点原理知牛顿迭代法对单根收敛.同理可知当*x 是0)(=x f 的重根时也收敛，则可分析出牛顿法不论对单根还是重根均是局部收敛的，只要初值足够靠近*x ，牛顿迭代序列均收敛于*x 。
四位二进制数可控加减法 上海交通大学电子技术实验大作业
《四位二进制数可控加减法》实验报告实验名称: 四位二进制数可控加减法姓名：学号：班级：目录一、实验方案 (3)二、设计思路................................................................................ 错误!未定义书签。
三、程序代码................................................................................ 错误!未定义书签。
四、调试问题 (6)五、心得感想 (7)一、实验方案1）基本功能实现两个四位二进制数的加减法运算,能够在led灯和数码管显示出结果。
2）清零功能利用一个微动开关，当微动开关按下时结果清零显示。
3）数码管显示将结果转换为七段显示器显示。
将运算结果输送到数码管中。
利用到人的视觉误差和短暂延时显示四位运算结果。
4）溢出问题若有溢出，则数码管显示“E”。
二、设计思路基本功能中分为连个模块，主模块用来运算加减法以及记录溢出和结果，子模块用来进行七段数码管的显示。
扩展功能中数码管显示要利用暂留现象，因此利用时钟clk来进行设计。
三、程序代码module show_sub(input [1:0]num,output reg [6:0] a_to_g );always @(*)case(num)2'b00: a_to_g=7'b1000000;2'b01: a_to_g=7'b1111001;2'b10: a_to_g=7'b1111111;2'b11: a_to_g=7'b0000110;default: a_to_g=7'b0000110;endcaseendmodulemodule show_top(input clk,clr,input wire [7:0] sw,input plus,sub,output wire [6:0] a_to_g,output reg [3:0] an,output reg [3:0] led );reg [15:0] clk_cnt;wire [1:0]s;reg [3:0] result; //运算结果reg [1:0] res;reg flag; //溢出标志wire [3:0] data1;wire [3:0] data2;assign data1=sw[7:4];assign data2=sw[3:0];assign s=clk_cnt[15:14];always @(posedge clk)beginclk_cnt=clk_cnt+1;endalways@(posedge plus or posedge sub or posedge clr)。
上海交大《计算方法》教学大纲
上海交通大学研究生（非数学专业）数学基础课程《计算方法》教学大纲(2007修改讨论稿)一．概况1.开课学院（系）和学科：理学院数学系计算数学教研室2.课程编码：3.课程名称：计算方法4.学时/学分：54学时/3学分5.预修课程：线性代数，高等数学，程序设计语言6.课程主干内容: 数值代数，数值逼近，非线性方程数值解，常微分方程数值解。
7.适应专业学科：全校的机、电、材、管理、生命和物理、力学诸大学科类，以及人文学科需要的专业。
8.教材/教学参考书：(1)李庆扬、王能超、易大义，数值分析（第4版），华中理工大学出版社， 2003(2)孙志忠，袁慰平，闻震初，数值分析，东南大学出版社，2002(3)J.Stoer and R. Bulirsch, Introduction to Numerical Analysis (secondedition), Springer-Verlag, Berlin-New York, 1993.(4)Atkinson K E，An Introduction to Numerical Analysis，John Wiley & Sons. 1989.二．课程的性质和任务本课程属于数值计算课程的基础部分。
数值计算课程是非数学类研究生数学公共基础课程，该组课程列入计算数学系列，目前按照“分级”的原则，设置《计算方法》（基础部分）、《微分方程数值方法》(扩展部分) 和《高等计算方法》（提高部分）三门课程。
本课程讨论用计算机求解数学问题的几类基本的数值方法及其相关的数学理论。
计算机是对近代科学研究、工程技术和人类社会生活影响最深远的高新技术之一，它对科学技术最深刻的改变，莫过于使科学计算平行于理论分析和实验研究，成为人类探索未知和进行大型工程设计的第三种方法和手段。
计算机的飞速发展正把计算的方法的创新、改进、提高推向人类科技活动的前沿。
人类现代计算能力的巨大更取决于计算方法的效率。
计算方法大作业
计算方法大作业班级XXXXXX 学号XXXXXXX 任课老师贺力平姓名XXX2013年12月实验一幂法与矩阵特征值1.幂法求主特征值思路幂法的主要思想就是对假设的任意初始列向量作用n次A矩阵（左乘A矩阵）后，初始向量就接近A矩阵的主特征值对应的特征向量。
由于左乘n次A矩阵有可能会造成计算量溢出，所以每次都对列向量作归一化处理。
2.程序代码function [ ] = mifa( A,v )%UNTITLED Summary of this function goes here % Detailed explanation goes hereA =[ 1 21 2 334 5 2 154 2 6 42 2 59 0];v=[1 1 1 1]';u(:,1)=v(:,1);fori=1:100v(:,i+1)=A*u(:,i);if abs(max(v(:,i+1))-max(v(:,i)))<10^-4breakendu(:,i+1)=v(:,i+1)/max(v(:,i+1));enddisp(u(:,i));disp(max(v(:,i)));k=u(:,i)'*A*u(:,i)/(u(:,i)'*u(:,i));disp(k);[x,c]=eig(A);disp(c);disp(x);disp(i);end3.结果比较和结论初始矩阵A =[ 1 21 2 334 5 2 154 2 6 42 2 59 0];初始向量v=[1 1 1 1]';幂法求得特征向量12.514714.914025.693539.6330归一化后特征向量0.31580.37630.64831.0000列向量最大值近似主特征值39.6330Rayleigh商求出主特征值39.6330用eig()函数算出的特征值和特征向量39.6331 0 0 00 -18.2401 + 7.4985i 0 00 0 -18.2401 - 7.4985i 00 0 0 8.8471-0.2450 -0.0471 + 0.0965i -0.0471 - 0.0965i -0.0574 -0.2919 0.1147 - 0.0939i 0.1147 + 0.0939i -0.1747 -0.5029 0.2862 - 0.1187i 0.2862 + 0.1187i 0.1535 -0.7758 -0.9330 -0.9330 0.9709 达到精度要求所需次数：17结论：可以看出初始列向量经过多次迭代后，用幂法求出的特征值和用eig()函数求出的A的特征值，满足计算精度在，并且特征向量也具有数乘关系。
2018年度交大计算机的第一次作业
1、多媒体计算机系统同一般计算机相比，特别需要利用计算机的数字化技术和______________ 。
选择一项：広a.通信技术D b.安全技术匕c.人机交互技术7匚d.存储技术2、计算机应用中，英文缩略语CAI所表示的计算机术语是_________ 。
* a.计算机辅助教学D b.计算机辅助设计匕c.计算机辅助制造匚d.计算机辅助工程3、当前气象预报已广泛采用数值预报方法，这主要涉及计算机应用中的____________ 。
匕a.数据处理和辅助设计□ b.科学计算与辅助设计-c.科学计算和数据处理d.科学计算和过程控制4、与二进制数110101110等值的十六进制数是屛b鼠标d.打印机/扫描仪6、假设给定一个十进制整数D,转换成对应的二进制整数B,那么就这两个数字的位数而言，B与D 相比， ________ 。
a. Db. Bc. DA d. B的位数大于等于的位数大于的位数大于的位数大于等于7、为解决某一特定的问题而设计的指令序列称为广a.b.语言文档c.系统卜d.程序丫8、对微机配置描述为中内存为 _________ 。
ra. 1.7G多媒体”，其b. 128MBI-c. 60GBa. 35ED b. 1AE fc. D70匚d. EA25、以下为计算机输出设备的是______________ 匸a.键盘d. 40X9、目前制造计算机所用的电子元件是广a.晶体管A b.集成电路电子管广c.超大规模集成电路10、有一类计算机应用领域称为 "CAD"，也就是b. 逢十六进 a. 计算机辅助教学c. b. 计算机辅助制造d. 逢八进15、八进制数631转成二进制数是c. 计算机辅助设计 计算机辅助管理d. a.10011101111、计算机内部采用二进制数进行运算、存储和控 制的主要原因是 广 b. 101011001a. 现/二进制只有0、1 两种状态，技术上容易实 c. 110011001 d.110100001二进制数数码少, b. 快 16、计算机的主机的组成部分是在计算机网络中传送速度 a. CPU 、外存储器、外部设备c. 二进制数数码少, 比十进制数容易读懂和记 广b. CPU 和存储器系统忆 二进制数数码少, d. c.主机箱、键盘、显示器存储起来占用的存储容量 小 12、在计算机中，所有信息的存储都采用 屛d. CPU 和内存储器17、下面哪一项不是计算机采用二进制的主要原因a. 卜六进制 二进制/b. a. 二进制只有0和1两个状态，技术上容易实c. 八进制 广bb. 二进制数的o 和1与逻辑代数的"真”和"假"十进制 d. 相吻合, 适合于计算机进行逻辑运算 13、计算机操作系统作为一个接口， 连接着广c. 二进制运算规则简单a. 用户与计算机7 金dd.二进制可与十进制直接进行算术运算b. 主机与外设 18、高级语言所编写的程序又称为源程序，此类程 序c. 系统软件与应用软件 a. 不大于100行的程序可以被机器直接执行 用户与软件d. A b.在更高级的大型计算机中能被机器直接执14、十进制数的运算法则是 a.逢二进 c. 不能被机器直接执行d.能被机器直接执行19、光盘是一种已广泛使用的外存储器，英文缩写CD-ROM旨的是_______ 。
上海交通大学管理学院《金融工程学》习题
一、大作业：本课程共包括3次大作业，旨在培养学生分析实际问题和解决实际问题能力。
要求学生自己实践与尝试，自己去调查、分析和计算，可以进行分组，进行学习小组交流、讨论，形成小组意见，课堂上安排小组代表作简要介绍，任课教师点评和总结。
1、设计“一个”新的金融产品。
2、计算一个具体的投资组合风险（例如VaR）以及解决风险的方法。
3、选择一个具体的金融产品定价（例如权证或者银行的理财产品）。
二、课后习题第１章金融工程概述1、请论述学习金融工程的三个基本目标，并举例说明。
2、根据已有的金融工程几个代表性定义，请阐述你对这几个定义的理解和看法。
3、请论述中国开展金融衍生产品交易的意义及其面临的问题。
第 2 章无套利定价原理1、假设市场的无风险借贷利率为 8 ％，另外存在两种风险证券 A 和 B ，其价格变化情况如图 2-11，不考虑交易成本。
图 2-11 两种风险证券的价格变化情况问题：（１）证券 B 的合理价格为多少呢？（２）如果 B 的市场价格为１１０元，是否存在套利机会？如果有，如何套利？（３）如果存在交易成本，例如，每次卖或买费用均为１元，结论又如何？2、假设无风险借贷半年利率 r ＝ 4 ％（单时期），两种资产的两时期价格变动情况如图2-12 ：图 2-12 两种资产的两时期价格变动情况问题：（１）利用动态组合复制定价技术给证券 B 定价；（２）如果证券 B 的市场价格为１００元，是否存在套利机会？如果有，如何构造套利策略？3 、试分析金融市场套利与商业贸易中的价差盈利的关系？为何金融市场中套利概念如此重要？第 3 章金融产品创新原理1 、如何设计一个更加合理的全流通方案？2 、如何设计一个金融新产品？第 4 章金融风险管理原理1 、金融风险是怎样产生的？如何从理论上解释金融风险？2 、怎样理解长期资本管理公司破产是一个由制度性缺陷、市场风险和流动性风险所造成的经典案例？3 、在例 4-1 中，当欧洲国家相关企业提出中国绍兴纺织企业向他们购买纺织设备，将终止使用美元支付的惯例，转为以欧元计价结算时，能否估计出１年之内因汇率波动产生的最大损失，若能是多少？4 、在例 4-1 中，能否找到一种套期保值方法，来减少思考题 3 估计出１年之内因汇率波动产生的最大损失。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
离步长改为 2.5mm*2.5mm。获得图像。
可见在 Bi 数较小的情况下，物体内的温度梯度较小，几乎可以近似地被认为是集总参数法 问题。 利用稳态条件下的计算， 以由外部热流流入的热量与由内部冷流流出的热量相同为计 算依据。 h1 S1 (T 1 T ) h 2 S 2 (T T 2) 。可得物体平均温度应当近似于 102 摄氏度。
T1 =141.7098 T2 =90.4344 T3 =84.0126 time =1.0257e+05 （2）取空间步长为 5*5mm，温度场矩阵的阶数为 20*10。Fo=6.25
7
T1 =144.0982 T2 =92.0991 T3 =85.0609 time =1.3781e+05 由数据可见，往往当使用更加大的微元体时，过程可以进行地更深入。这不难理解,首先， Fo 数是代表非稳态进程程度的无量纲时间数， Fo 更大， 过程进行地越深入。 另一方面来讲， 因为更长的步长，意味着更加长的时间步数。在相同的收敛判断条件下，由于步长更长的运 算单步的温度变化更大，所以进行的程度更深。
2
3
T1 =133.6058 T2 =80.6383 T3 =74.5586 time =5.3573e+04 设导热系数为 40w/m*K,内部空洞的传热系数为 3W/ m 2 *K， （2）使导热系数大幅增大，
2 外部传热系数为 5W/ m *K。计算可得 Bi1=0.00125,Bi2=0.00075。为了缩短计算时间，将距
3.时间步长对结果的影响
（1）.首先取时间步长dt=x^2/(a*4)=1.875，Fo=0.25，
T1 =133.6058 T2 =80.6383 T3 =74.5586 time =5.3573e+004得到标准数据
8
(2).取时间步长dt=2，Fo=0.267，发散
(3).取时间步长dt=1.7，Fo=0.227，
其中 i , j ∈ Ζ，Fo = a ∗ dt x 2 ，Bi = hx λ，a = λ ρc，dt = x 2 4a
三、编程思路
把 50*100 的节点的温度当做一个 50*100 的温度矩阵。因为在非稳态显示格式中，大多 数的节点的热平衡方程都是相同的。 所以把这个矩阵的全部内容全部往四个方向移动一行或 一列，获得四个新的矩阵。再根据热平衡方程，用这四个新的矩阵和原矩阵进行矩阵加减运 算，就可以得到一个时间步长后各个节点的温度矩阵。 在有边界的地方（M，N，P，Q）依然按照类似方法制造四个矩阵，但是有边界条件的 方向使用边界条件给定的特殊值来创造矩阵。 然后同样利用显式方程来计算一个时间步长后 的温度矩阵。 得到一般点和所有点的温度之后， 再把各个矩阵中有用的部分拼接起来就得到了完整的一个 步长之后的温度矩阵。 利用循环语句反复进行上述步骤，直到判定收敛。 程序代码见本报告附录。
Bi h dx /
目标是借助 matlab 用数值分析的方法求得整块板达到最终稳态的过程中， 这块板的温度 随时间的变化，以及其中三个代表性点的温度随时间的变化，并通过调整不同的时间步长， 空间步长来探索其对结果的影响，从而掌握数值方法的运用，进而验证传热学的理论知识， 通过实践及思考获取新知识。 由图可得，图形左右对称，即轴线位置属于绝热边界，所以我们可以把板沿轴线左右一 分为二，只要计算一半的板即可，这使需要的计算的范围缩小了一半，降低了工作的难度，
T1 =132.6748 T2 =79.4658 T3 =73.3294 time =5.0849e+004 当稳定性条件都满足才能保证整个方程组稳定，本题情况下需要保证1-4Fo≥0，当1-4Fo ＜0的时候结果不稳定。满足1-4Fo≥0时，当x一定，dt越小，Fo越小，越能满足稳定性 条件。但是相应的计算工作量增大。
2 2
6
可见结果与上一种情况完全相反。 T1 =31.7660 T2 =23.6478 T3 =22.8033 time =15540 因为从初态到稳态的温度变化很小，所以花费的时间也最小。
2．空间步长对结果的影响
为了使结果具有一般性，使用在自然空间较为常见的导热和传热系数。设导热系数为 40w/m*K,内部空洞的传热系数为 30W/ m 2 *K，外部传热系数为 50W/ m 2 *K，外部流体温度 为 200 摄氏度，内部流体与对象初始温度为 20 摄氏度。 为了便于阅读，将标准计算的数据列于此。 T1 =133.6058 T2 =80.6383 T3 =74.5586 time =5.3573e+04 （1）首先取空间步长为 2.5*2.5mm，温度场矩阵的阶数为 40*20。Fo=1.5625e-04
2 2
5
易于理解，在内部流体的传热系数较小时，整个传热过程被外部流体所主导。内部流体需要 更高的边界层温度梯度，来吸收外部流体所带来的热量。 T1 =172.1985 T2 =151.5862 T3 =148.6036 time = 1.0529e+05 因为初始条件与稳态之间的差别较大，导致消耗了很多时间。 设导热系数为 40w/m*K,内部空洞的传热系数为 （5）使外部流体的传热系数大幅减小， 30W/ m *K，外部传热系数为 5W/ m *K。计算可得 Bi1=0.00125,Bi2=0.0075,。获得图像。
1.导热与传热参数对结果的影响
为了使数据有一般性，选取 1mm*1mm 的控制体，则温度场矩阵的阶数为 100*50。外 部流体温度为 200 摄氏度，内部流体与对象初始温度为 20 摄氏度。 设导热系数为 40w/m*K,内部空洞的传热系数为 30W/ m 2 （1）选取一般性的常见参数， *K，外部传热系数为 50W/ m *K。计算可得 Bi1=0.0125,Bi2=0.0075,Fo=0.000025。获得图像。 （这一次运算条件作为标准计算，将会被多次提及）
4.与稳态解的比较
稳态过程的计算本质上与非稳态有相似之处，在此不加说明，使用代表稳态运算的 m 文件可以得到与第一次标准状态相同的条件下的解。 T1 =110.8181 T2 =79.4591 T3 =74.8974 显然与标准运算有比较大的出入，原因不明。
9
附录
1．data
%设定节点长度，假设整个图形是一个 1*0.5 米的区域，x 为单个微元的长，y 为宽。 %x 与 y 取 0.01m，0.025m，0.05m 等，可以使长度方向微元数为 4 的倍数，宽度方向为 2 的倍 数。 t0=20 tf=200 h1 = 50 h2 = 30 pc = 30e5 x = 0.1 y = 0.1 gama = 40 %计算得到的常量 a = gama/pc dt = x^2/(a*4) Fo = a*dt/x^2 Bi1 = h1*x/gama Bi2 = h2*x/gama %辅助量，X为长度方向的微元数，Y为宽度方向上的微元数 X = 1/x Y = 0.5/y A =ones(Y,数值解相近。边界层的梯度非常大，可见在这种情况下，总的传热量最大。 T1 =108.2997 T2 =100.5442 T3 =99.3991 time =3.8500e+05 因为传热系数小，系统吸收热量的速度很慢，所以从初态到稳态花费的时间很多。 （3）使用非常小的导热系数，设导热系数为 4w/m*K,内部空洞的传热系数为 30W/ m 2 *K，
计算方法期末大作业
一、问题分析
已知量 t 0 = 20℃ t f = 200℃ h1 = 50W m2 ∗ K
h2 5W / m2 * K
ρc = 30 ∗ 105 dx = 0.01m λ = 60W m ∗ K 导出量
a / c
dt (dx) 2 /(a 4)
Fo a dt /(dx) 2
10
A3=[B D]; B=A(1:49,1:100); C=A(1:1,1:100); A4=[C;B]; G=(A1+A2+A3+A4)/4; %然后计算有对流换热的边界的点，首先计算M边 M=A(1:50,100:100); M1=A(1:50,99:99); C=A(2:50,100:100); z=A(50,100); M2=[C;z]; C=A(1:49,100:100); z=A(1,100); M4=[z;C]; C=ones(50,1); M3=tf*C; D=gama*(M1+M2+M4)/(3*gama+h1*x)+(h1*x*M3)/(3*gama+h1*x); M=D; %然后计算N边 N=A(26:50,76:76); N3=A(26:50,77:77); N4=A(25:49,76:76); C=A(27:50,76:76); z=A(50,76); N2=[C;z]; C=ones(25,1); N1=t0*C; D=gama*(N3+N2+N4)/(3*gama+h2*x)+(h2*x*N1)/(3*gama+h2*x); N=D; %计算P边 P=A(25:25,26:75); P4=A(24:24,26:75); P3=A(25:25,27:76); P1=A(25:25,25:74); C=ones(1,50); P2=t0*C; D=gama*(P1+P3+P4)/(3*gama+h2*x)+(h2*x*P2)/(3*gama+h2*x); P=D; %计算Q边 Q=A(26:50,25:25); Q1=A(26:50,24:24); Q4=A(25:49,25:25); C=A(27:50,25:25); z=A(25,50);
1
也为后续的编程工作节省了计算的时间：
要运用数值分析的方法， 首先我们要将这块板的这块区域划分成许多互不重叠的单元体。 每个单元体的物理量用某一点的物理量来代替， 然后利用每个节点的物理量都与它周围的节 点的物理量有一定关系，将物理现象描述成一系列代数方程组，再用计算机求解。首先，由 于此题属于非稳态导热，选用 x-y-t 的三维坐标系。然后，将 50*100 的区域划分为 50*100 的节点，并列出相应节点的节点方程。最后，用矩阵法直接联解代数方程组。