[讨论]

偏微分方程式数值解之计算机算法详解

2020-8-29 14:31:01 2020

0 偏微分方程式数值解之计算机算法详解 ---Fortran语言数值方法编程实现内容与附图内容一致，符合国际标准，附图待续。 ***************************************************************** 本章节偏微分方程（计算机数值方法语言Fortran解答） Pg.lp343 导论框架1 ---- ----- 偏微分方程常常发生在实践中。有些例子如下：位移y，一根紧紧延伸的绳子经历横向振动，给出。。。？（方程式）？。。。，t为时间， x为从一端的距离，且c为一个常数。一均匀的线束经历横向振动的位移类似给出为。。。？（方程式ii）？。。。两维度稳态热流方程式是。。。？（方程式iii）？。。。.这个方程式给出了？西塔O/点的温度，在一个平的盘子上其坐标为？（x，y），它是隔热的没有热量流进流出盘子，并假设过去了足够的时间没有进一步随时间变化。若温度也随时间变化，那么方程式成为。。。？（方程式iv）？。。。我们就此不考虑这些方程式的推导，但若您感兴趣您将在我们的“工程师和科学家的数学”第二卷中找到一定数量的偏微分方程的推导。 ------------------------------------------ 框架 2 ------ 举个例子，在上一个框架中给出的第一个方程式，很容易见到？y=f（x+ct）和？y=g（x-ct）是解析解答，这里f和g代表的是任意函数。进一步讲，？y=f（x+ct）+g（x-ct）也是一个解。验证是如此通过微分和替代。 ************************** 2A -- 。。。？（偏微分方程式组v）？。。。所以结果。 ------------------------------------------ 框架 3 ------ 虽然好像若存在这样一个解析解则我们所有的麻烦都没有了，情况绝不是这样的。当我们试图满足初始条件时困难出现了，如所成立的相当不适合。这就是为什么使用了像变量分离的特殊方法。但是甚至这些特殊方法不总能够帮的到，而在这样的情况下尝试了数值技巧。在此程序中我们将给您显示仅一种方法处理这样的问题的简单例子。 ------------------------------------------ ************************** 偏微分方程 Pg.lp344 框架 4 ------ 拉普拉斯方程式。。。？（偏微分方程式）？。。。 --------------------------------------------- 假设一个方形的盘子的表面隔热没有热量流进流出它们。让盘子的一边保持在100摄氏度而其它三边在零度，这些温度精确到接近十分之一度。经过了足够的时间之后，盘子上的任何点的温度将独立于时间。现在这个问题可以通过变量分离的方法解答，但是取而代之的将是足以用于方便地演示数值技巧。 ---------------------------------------------------------- 框架 5 ------ 开始做的第一步是把方形分成数个较小的方块。除非有理由不这样做，边缘好像首先划分成两块，然后四，接着八，等等。那么产生了一个网格，而这种方法包括计算小方块的角上的温度。示意图仅仅显示了网格的一部分。让每一小方格的边是h。那么取网格点P的坐标为？（x，y），Q，R，S和T的坐标相应的坐标是（x+h,y）,(x,y+h),(x-h,y)和（x，y-h）。。。。？（SPQT网格示意图）？。。。 --------------------------------------- 框架 6 ------ 现在采用在前回程序框架52使用的过程。这里在点？（xn，yn）的？d2y/dx2由？（算式）？代替。您可以写出在？P点的？d2O//dx2和?d20//dy2对应的公式吗？写下P点?O/的值为？O/(x,y),等等。 ********************************* 6A -- 。。。？{算式12}？。。。 ---------------------------------------------------------- 框架 7 ------ 代替这些到方程式？（方程式）？给出，估计，。。。？（算式）？。。。（7.1），例如，。。。？（算式）（7.2）。这说明在任何点P的值近似等于在？Q，R，S和T四个点的平均值。 ———————————————————————————— ************************** 偏微分方程 Pg.lp345 框架 8 ------ 。。。？（0到100度，ABDEFGHJK示意图）？。。。现在假设原始方块首先一划为四，然后为16，小的方块。让在小方块角的温度如图所示。下一步是找到所示的九个温度的初始估值。起初找到？西塔O/22的估值。这样做是通过把P，Q，R，S，和T在A，B，D，E，和F.那么？西塔O/22的初始估值是在B，D，E，和F的平均值，例如，1/4（100+0+0+0）=25。 --------------------------------------------------------- 框架 9 ------ 下一个估算的值是？西塔O/11. 而？西塔O/12和？O/21还不知道，（7.2）不能直接使用。代之，一个斜的方块（显示点划线）画在里面，且？西塔O/11的值取四个值的平均值，在这个方形边的中点，例如，在A，D，G，和E的平均值。这样得出初始估算的？西塔O/为1/4（25+0+0+0）=~6.2. 至1小数位，什么将是？西塔0/13,31和33的初始估算？假设在H和J两者的温度是？50摄氏度。 ************************** 9A -- 。。。？（？西塔O/13,31,33=数值序列）？。。。 ---------------------------------------------------------------- 框架 10 ------- 直接使用（7.2）现在可以得到？西塔O/12,21,23,和32的估值。这些的估算值将是什么？ ****************************** 10A --- 。。。？（？西塔O/12,21,23和32=数值序列）？。。。 ---------------------------------------------------------------- *************************************************** 偏微分方程 Pg.lp346 框架 11 ------- 我们现在首先有估算在网格点的所有的九个温度。一种高斯-塞戴尔重复技巧现在用于改进这些值。因为这是可能的，由于您将能过很快验证，关系系数矩阵是斜支配的。用（7.2）继续写下所有九个？西塔/的它们附近的一般表达式，所以。。。？（？西塔/1112,22）?。。。其它？西塔O/的对应方程式将是什么？ **************************** 11A --- 。。。？（？西塔/13...32=算式阵列）？。。。在每种情况下右边的四项写为逆时针次序，从找到的温度点。这样做组成了一种系统的方式列出方程式。 ------------------------------------------------------- 框架 12 ------- 通常这些九个方程式现在会依次使用来改进？西塔/值。在这个特别的例子中，关于？EB明显有对称，所以有？西塔31=11，？32=12和33=13。利用这一对称减少九个方程式。。。？(西塔11...23, (12.1),(12.2)算式阵列)?。。。 -------------------------------------------------------- 框架 13 ------- 在（12.1）使用初始估值？西塔/12和21得出第二次估算？西塔O/11为（18.8+9.4）/4=7.0. 使用这个值一起与？西塔13,22得出12的第二次估值为（43.8+7.0+25）/4=19.0. 其它四个？西塔O/的第二次估值将是什么？在每种情况下使用现有的最近的值来计算。 ****************************************** 13A —- 。。。？（？西塔O/13,23=数值算式）？。。。 -------------------------------------------------------------------- 框架 14 ------- 现在的过程按必要重复许多次直到估算上没有发生进一步的更改。做工（手工进行时）最好列在一张表内： ************************** 偏微分方程 Pg.lp347 框架 14（续） ------------ ？西塔O/11,12,13,21,22,23 ------------------------- 。。。？（数值表格）？。。。得出这些值到小数一位，如它开始的那样，引用的原始温度纠正到十分之一度。我们到达阶段时停止计算，计算的最后六个数值，在一对一的基础上等于前回的六个。这意味着若再进行计算值将不会进一步改变。当使用一台计算机时，您可能会指引它测试计算到每列最后的收敛。在那样的情况下，输出会包括额外的读数25.0和52.7，在没有完的列的尾部。 ----------------------------- 框架 15 ------- 通过把目前的小方格化成四个，原始的方形现在会划分为更小的方形。那么将有64个小的矩形。若做完的话，两件事将会发生。这些的之一是温度将从49个内部的点获得而不是仅仅9个。第二是您会期望结果的更大的精度，一般在数值工作中， h值越小，结果越准确。由于减少小方格的大小没有什么新鲜事儿，我们将不最求结果的计算。 ----------------------------- 框架 16 ------- ？西塔O/11...示意图？。。。？（0到100度）？。。。示意图展示了刚刚做的类似的问题，但是这次是一个矩形区域划分为数个小方格。对于这种设立的那些框架11和11A中类似的方程式将是什么？ **************** 16A --- 。。。？（西塔O/11...算式阵列）？。。。 ----------------------------- ************************** 偏微分方程 Pg.lp348 框架 17 ------- 在前回的例子中可能获得重复的合理准确的初始值，且若可以这样做的话，最好就这样做。这里就如此简单，对于一个开始，矩形的中点不在网格的点上。可是，按常理您可以写下一组值，并从这一组开始重复。现在明显，顶行的？西塔O/值明显将大于底部的那些。类似地您将期望在示意图右边的那些将大于左边的那些。所以开始的一组可能的数据取为。。。？（？西塔/11,...24=数值阵列）？。。。使用这些数值和？16A中的方程式，找到?西塔O/的值，取最近的整数。 *********************************** 17A --- 。。。？（西塔O/11...24数值表格）？。。。 ------------- ----------------------------------------------------- 框架 18 ------- 如Pg.349所示一张流程示意图，拉普拉斯方程式的解答给一个矩形网格，横向n个间隔，及m个向下。 ----------------------------------------------------- 框架 19 ------- 这类例子算法的另一种做法是放大示意图并在图中插入找到的连续的值。这得出一个更为自动的过程，当手工做题时需要较少的思考，但不适合计算机使用。示意图如Pg.350所示，演示了框架16中例子的理念。数据计算实际上于17A中一样次序，但是代替了在这个表中查找以前值插入相关的方程式，找到的每一个值取到最右边值的平均值，到最左边，和最下边。所以例如，双星的64是四个单星数的平均值，例如， 63,100,60和32。 ------------------------------------------------------ 框架 20 ------- 方程式？（方程式=0）？已知为两个维度的拉普拉斯方程式。本程序中将要处理的其它方程式是。。。？（偏微分方程式i，ii）。 ************************** 偏微分方程 Pg. 349 框架18流程图。开始读入读入边界值由于计算中不需要“角”值，它们可以缺省。在内部的网格点初始化值所有值从零开始可以更容易试图获得首次估算。收敛会变得较为长一点，但是如使用计算机差别会忽略掉。。。。？（算式）？。。。收敛发生在，？不同于前一个值少？。是否停止（解答。。。？偏方程？。。。参考（4,5）。） ************************** 偏微分方程 Pg.lp350 。。。？（框架19.示意图）？。。。框架20（续） ------ 基本上我们给这些的每一个使用的方法是将类似于给拉普拉斯方程式。可是有一些点要注意关于这些其它方程式。 --------------------------------------------------- 框架 21 ------- 方程式。。。？（偏微分方程式）？。。。 ----------------------------------------- 当解拉普拉斯方程时，？h用于两者的增量和y。这是合理的，由于x和y均为距离，且相等的增量发生当x，y图划分为数个小方格。现在？（方程式）？是一根振动的琴弦的方程式，且这里x代表一段距离，而t代表时间。所以没有 ************************** 偏微分方程 Pg.lp351 框架21（续） ------- 特别的原因给这两个变量取相同的增量。在这这种情况下让增量在x和t是h和k对应。现在，对于这个微分方程式，对应于（7.1）的方程式将是什么？ ****************************** 21A --- 。。。？（算式）？。。。（21A.1） -------------------------- 框架22 ------ 这个方程式取一个特别简单的形式，如果？h和？k之间取一定的关系。记得x和t是独立的变量所以我们可以取适合我们自己的h和k。您能够建议最佳关系h和k以使（21A.1）成为尽可能简单？这一简单格式是什么？ ******************************* 22A --- 如果？（算式），例如，如果？h=ck，（21A.1）成为。。。？（算式）？。。。（22A.1）或者这样的另一种版式，其中h或者k不出现。（您假设h和k两者都是正的.） ---------------------------- 框架23 ---- 所以，。。。？（算式）？。。。（23.1）且这是我们将要使用的这个方程式的格式。要演示，让我们取方程式。。。？（算式）？。。。并找到一个解答其满足条件。。。?（算式）i）ii，iii，iv））？。。。对于 0《x《2. 现在已经开始微分方程是一根振动的琴弦。运动的信息及给出的条件它是如何开始的给了您什么信息？ ********************* 23A --- 点？x=0和？x=2永久是静止在零位移。运动开始通过绑定琴弦为半个正弦波并从那个它静止的位置释放。 ------------------------------ ************************** 偏微分方程 Pg.lp352 框架 24 ---- 现在假设x的增量是0.2。t的增量您取多少？ ******************** 24A --- 0.1，从h=ck有h=0.2和c=2 ------------------------------- 框架 25 ------- x的值局限于0和2包含在内之间。t的值将无限地增加。这意味着，思考到时间，仅仅显示表格的开始。然后可以按必要扩展。一些信息可以马上插入表格，如下：。。。？（t，x数值表格）？。。。所有表格总体内的值乘以10-6. 在行t=0的值是那些函数1/100sin1/2πx给0（0.2）2。您将注意关于x=1是对称的，且这在表中继续。 **************** 25A --- 。。。？（t=0方程式）？。。。 ----------------------------------- 框架 26 ------- 现在得出下示点状表的一部分。。。？（点阵示意图）？。。。方程式（23.1）给出一个公式给方形点代表圆点，如位置所示，挑出的点相对应， ************************** 偏微分方程 Pg.lp353 框架 26 （续） ------------- 。。。？（算式阵列）？。。。所以当填上数条线时，下一个找到通过移动方形的点沿着现在做的一条线。如果我们试图以这个步骤开始表格，您可能发现多么的难发生了？ ***************** 26A --- 做出行t=0.1的输入项，在每次情形下，最上面的点离开点。 ------------------------------ 框架 27 ------- 找到这个难的行的输入，我们必须包含目前没有使用的初始条件，例如，？（dy/dt）t=0=0. 这是通过下列方式完成：在我们的实际问题中，运动以独特的方式在时间t=0开始。但是，若不从这个时间开始，让我们假设发生在t=0之前的方法是，在t=0，琴弦实际上是按问题讲的它应该的。当？t=0时，（22A.1）成为。。。？（算式）？。。。（27.1）且，虽然我们不知道？y（x，-k）的值，假设其与我们的“发明”一致。 ------------------------------------------------ 框架 28 ------- 现在，？dy（x，t）/dt的一个近似的表达式是？(算式)，所以。。。？（方程式）其中，当？t=0时，给出。。。？（方程式）？。。。. 这个结果现在包含在（27.1）. 当这个完成时，您得到了y（x，k）啥？ ******************************* 28A --- 。。。？（方程式i）？。。。所以。。。？（方程式ii）？。。。（28A.1） ------------------------------------------------ 框架 29 ------- 这个方程式现在能够使我们填入难的行。在特别做的例子中，？（方程式）和这里（28A.1）成为。。。？（算式）.?。。。？（方程式i），所以这里（28A.1）成为。。。？（方程式ii）？。。。 ************************** 偏微分方程 Pg.lp354 框架 29 （续） ------- 所以当x=0.2时，在行t=0.1的输入，例如，y（0.2,0.1），是1/2（0+5878）且当x=0.4时是1/2（3090+8090）. ************ 29A --- 。。。？（数值列表）？。。。 ------------------------------ 框架 30 ------- 这之后，（22A.1）直接用以给出输入，为了下列行。所以，当t=0.2，x=0.2，输入得出0+5590-3090=2500，且当t=0.2，x=0.4时，对应的输入是2939+7694-5878=4755。现在完成输入行给t=0.2. *************************** 30A --- 。。。？（数值列表）？。。。 ---------------------------------- 框架 31 ------- 现在可以填入表内，当t=1.0时，成为。。。？（t，x数值表格）？。。。不要忘记每个输入乘以10-6。另外，由于在框架25上提到的对称，仅仅显示了x=0和x=1之间表的部分。 ---------------------------------- 框架 32 ------- 刚刚解决的问题可以解析地处理通过分离变量，且数值技巧不是必要的。解析解是？y=1/100sin1/2πxcosπt及，若您希望，您可以插入上表给出的x和t的任何对子值到这个方程式，并且检查数值解答是多么的准确。 ---------------------------------- **************************** 偏微分方程 Pg.lp355 框架 33 ------- 下列问题您会发现有麻烦，若您试图采用变量分离的技巧解答它：一根延伸的琴弦处于静止的位置如图示。当震荡时，运动满足方程式？（二级偏微分方程式）？。运动传递到琴弦通过给出一个强迫的震荡？y=1/50cos2πt到端点A，其它的端点保持固定。问题是在后续的时间找到其它点的位移。。。。？（01xy1/50A示意图）？。。。若您熟悉变量分离的方法，以这种方法试一下看看发生了什么。否者直接进到下一个框架。 ****************** 33A --- 分离变量导致。。。？（算式）？。。。若 y=0 当 x=0，A=0。当？dy/dt=0 当t=0，？B1=0。所以。。。？（y算式）？。。。当x=1，。。。？（y算式）？。。。，所以？k=1/2π，所以？C=1/50。。。。？（y=1/50算式）？。。。当t=0，？y=1/50x 及不能满足这一条件。 --------------------------------------- 框架 34 ------- 若x取在间隔01，t将取在什么间隔？ ******************** 34A --- 0.025 ---------------------------------------- 框架 35 ------- 构建一个类似于那个在框架25的表格，填到顶行和第一及最后一列。取t值从0到0.5，并做到5位小数。 ******************* ************************** 偏微分方程 Pg.lp356 35A --- 。。。？（tx数值表）？。。。所有表的主体的输入乘以10-5。最后的一列通过使用函数得到？1/50cos2πt。 ------------------------------ 框架 36 ------- 现在完成表格，如前回做完相同的线。 ****************** 36A --- 表如357页所示。 ------------------------------ 框架 37 ------- 如果您检查表中的数据您刚刚完成的，您将注意到一或者二个有趣的点。首先，之前一些时间运动的影响传递到端点A沿琴弦到达另一点。第二琴弦上的一些点的位移在一定的时间上，超出了端点A的最大位移。第三点的最大位移随时间变化，所以给出了一种波的效应。 ------------------------------ ************************** 偏微分方程 Pg.lp357 。。。？（框架36A. tx表）？。。。框架 38 ------- 解答？（方程式）？的一种流动图有？dy/dt=0在？t=0，y固定在两端对于所有的t值，如框架23，在358页所示。 ------------------------------------ 框架 39 ------- 方程式？（方程式） ----------------- 这是控制一根杆子上的热流的方程式，其侧边是隔热的。？西塔O/温度是距离和时间的一个函数，且再讲h和k的增量将对应地取在x和t。对于这个解答的数值解，使用了如前回的同样类型的表达形式，给二级导数和？{算式}？给一级。在这种情形下方程式对应于（7.1）将是什么？ *************************** ******************************** 偏微分方程 Pg.lp358 框架 38 ------- 。。。？（流程图）？。。。开始读入读边界值计算输出值输出值是否停 m=t的间隔数 n=x间隔数在t=0时，边界值可以由一个公式指定（例如），其中必须计算。。。。？（算式12）？。。。这个程序流程图显示在编程中使用双字符如何是没有必要的，若每一行的值在计算时打印出来了。在任一时间仅仅需要保存三行的值。在计算机Fortran语言中这些可以标记为Y012（JJJ）。对这张程序流程图略作修改会适合框架33情况的类型，这里一个端头不固定。（解答方程式。。。？（偏微分）？的程序可以在参考（245）中找到）。 ******************************** 偏微分方程 Pg.lp359 39A --- 。。。？（算式）？。。。 ----------------------------------- 框架 40 ------- 重新排列后，这可重新列为。。。？（算式1）？或者。。。？（算式2）？。。。（40.1）若?kc2/h2由？阿尔法a标记。现在可以理论上显示？阿尔法a将是正的，应该小于或者等于1/2，否者解答中可能建立大的误差。这里我们将不投入到这一理论，但是您若感兴趣，您将会在章回8中找到处理，“微分方程的数值解”，由W.E.米尔尼（威利）著。进一步可以显示，从精度的观点上，给？阿尔法r取的最佳值是1/6。当这项完成时，（40.1）成为当这项做完后，（40.1）成为。。。？（算式）？。。。（40.2）这里？h和？k的连接关系式为？（算式）？. ------------------------------------- 框架 41 ------- 看一下这是如何实现的，假设一根隔热的杆子长度1米，一端头保持在100摄氏度，且另一端在零度。容许足够的时间流逝以使稳态条件流行。热端的温度突然减到零。假设？c2的值是1/60，温度将在后续时间找到在10米的间隔沿着杆子的长度。为了使用这一方程式（40.2），什么将是k要求的值？（记得所有的长度将由米数表达。） **************************** 41A 。。。？（k=算式）？。。。 ----------------------------------------------- 框架 42 ------- 当稳态条件流行时，沿杆子10厘米的间隔处的温度将是0（10）100. （温度将从一端向另一端均匀增加。）假设在热端头的温度的降低发生在时间？t=0。所以在立刻？t=0之前，？西塔O/在这一端是100，而紧接着，它是零。可是这一改变占有了一个有限的时间，所以合理地取这一端头的温度为50摄氏度在？t=0，且紧接着零。所以我们的初始值表格将是： **************************** 偏微分方程 Pg.lp360 框架 42（续） ------- 。。。？（t，x表格）？。。。 ------------------------------ 框架 43 ------- 我们现在可以开始在行？t=0.1填充值。当？x=0.1,使用（40.2）得出？西塔O/（算式）=10.类似地，当？x=0.2，？西塔O/=(算式)。这一行将要填充的剩余的值是啥？（按必要做到1小数位）。 **************************** 43A --- 。。。？（数值列表）？。。。 -------------------------------------------- 框架 44 ------- .................... 。。。？（t，x数据表）？。。。 ******************************* 偏微分方程 Pg.lp361 框架 44 （续） ------- 若您检查前页表内的数据，您将见到它们显示了您会期望的非常清晰准确的影响。 ----------------------------------------------------------------------- 框架 45 ------- 。。。？（方程式）？。。。的解答的一张流程图给一维度的流动系统有常端头条件显示如362页。？x和？t是这样的？0< ----------------------------------------------------------------------- 框架 46 ------- 对于上个问题的稍微的修改示范按条件改变如下所示：让端头？x=0隔热并假设热量从杆子的另一端头丢失当流体流过时。正如这样，流体将当然变得发热。让杆子的初始温度是1000摄氏度，且流体流动的温度是50摄氏度，在与杆子接触之前。这里与上一个问题的差别是端头冷却，其温度没有突然冷却至最终值，而温度慢慢地发生变化。为了包含这种慢慢的改变，有必要引入温度梯度？dO//dx?在这个端头。 ----------------------------------------------------------------------- 框架 47 ------- 现在可以显示，在这些情况下，从一个表面丢失的热量由两个不同的表达式之一给出： i）？HA（？西塔O/-O/c）这里？H是传热系数， A表面的面积，？西塔O/表面的温度，及？西塔O/c流体的温度接近表面， ii）？-KAdO/dx这里？K是杆子的到热率。把这两个表达式相等，。。。？（方程式）？。。。（47.1）对于？x=L，L是杆子的长度，这里是1米。 ------------------------------------------------------------------------ 框架 48 ------- 为了使用这个方程式，采用了一种有些类似于框架27内的灵巧方法。我们想象杆子增加了?h(0.1米这里)的长度，且温度在？x=1.1以这样的方式变化，即在实际端头温度总是在实际问题的条件应该的之下。在实际端子的温度现在使用（40.2）来计算。使用？L和？h而不是1和0.1（为了得到一个更一般的结果），（40.2）得出，对于？x=L 。。。？（算式）？。。。（48.1）且，在每一阶段，未知量在右手边是？西塔O/(L+h,t). ----------------------------------------------------------------------- ************************ 偏微分方程 Pg.lp362 框架45流程图。。。。？（程序流程图）？。。。开始读步长步数本文中，？阿尔法a固定后找到？k 程序中结合检查值可接受 y在t=0的分布这些是常数和条件。计算输出停止是否注意在编程中，不必要使用双字符。一次仅需要存储两行的？y值。在计算机Fortran语言中这些标记为YO（J）和Y1（J）。 ************************ 偏微分方程 Pg.lp363 框架 49 ------- 通过使用（47.1）寻找未知量。？dO/dt,对于？x=L和时间？t，可以由？（算式）？逼近，所以？（47.1）成为，对于？x=L，。。。？（方程式）？。。。（49.1）包含这个方程式到（48.1）以致结果是？西塔O/(L,t+k)的一个方程式，不包含？西塔O/(L+h,t). ******************************* 49A --- 。。。？（算式）？。。。（49A.1） --------------------------------------------------- 框架 50 ------- 现在我们的？h值是0.1，且？k也是（从？41A），并取按必要？H和？K的典型合理的值，3，000W/m2摄氏度和200W/m摄氏度对应，（49A.1）成为，有？L=1和？西塔O/c=50, 。。。？（算式）？。。。（50.1） (40.2)现在用的值是？x=0.1(0.1)0.9，且这最后的方程式对于？x=1。 --------------------------------------------------- 框架 51 ------- 我们还有问题寻找连续的？西塔O/值，当？x=0时。这样做的一个过程类似于刚刚给？x=L采用的。可是结果，有点简单，由于在端头，？x=0，？dO/dx=0,如杆子这里隔热。通过适当地修改（48.1）和（49.1）您可以做出公式，对于？O/(0,t+0.1)? ***************************** 51A --- （48.1）将被下式替代。。。？（算式）？。。。（51A.1）（49.1）将被下式替代。。。？（算式），给出？（算式）（51A.1）那么成为。。。？（算式）？。。。例如。。。？（算式）？。。。 ------------------------------------------------------ ********************************* 偏微分方程 Pg. lp364 框架 52 ------- 三个方程式为。。。？（算式123）？。。。现在及时用于连续的每一阶段。对于？x从？西塔o/(x,0)=1000开始，给这个新的问题对应于框架44做出这张表。直到t=1.0连续解答，给出?o/值到最近的整数。 **************** 52A --- 。。。？表x t？。。。 ------------------------------ 框架53 ------ 计算温度分布的一个流程图，以连续的时间，在一个一维度的热流系统一端？（x=0）绝热及一个热增益（或者丢失）在另一端？（x=L），如lp365页所示。微分方程组为：。。。？（偏123）？ ------------------------------ 框架54 ------ 方法中的误差 ------------ 在这个程序中处理的问题的解答，不建议深入到误差的细节分析。这些会由于两种原因。第一是由于代替偏导的表达式例如。。。？（算式），？。。。它们本身仅仅为近似的。这样的效果是替换偏微分方程由已知所谓的一个差分方程。 ********************************** 偏微分方程 Pg. lp365 框架53的流程图。。。？（示意图）？。。。开始读取计算步长计算m和n 这是在t=0时系统中的温度分布输出是否停止（解答。。。？（偏微分方程）？。。。的程序可以在参考（2,4和5）中找到。） ************************************** 偏微分方程 Pg. lp366 框架 54（续） ------- 一般地讲，您不要期望差分方程的解答具有它替代的原始的偏微分方程的相同的答案。而不寻常在情形方程式。。。？（偏微分方程）？。。。，若？h和？k满足关系式？h=ck，两个答案是相同的。所以简化差分方程的同时，这也是考虑精度取的最佳值。验证?y=f(x+ct)+g(x-xt)满足？y(x,t+k)=y(x-h,t)+y(x+h,t)-y(x,t-k) 若？h=ck. ********************* 54A --- 在差分方程式中替代左手边=。。。？算式？。。。右手边=。。。？算式？。。。=左手边 ---------------------------- 框架 55 ------- 现在您已经在2A中验证了这个解答满足。。。？（偏微分方程）？及在条件？h=ck下，差分方程和偏微分方程均有相同的解答。 ---------------------------- 框架 56 ------- 误差的另一个源是由于算术取整。甚至这类误差可以导致严重的麻烦，情形是在方程式。。。？（偏微分方程）？，如果取？阿尔法a大于1/2。 ---------------------------- 框架 57 ------- 在这个Fortran程序中，我们已经集中精力给你们展示了一种基本方法处理一定的偏微分方程。另外还考虑了其它一些方法和升华在基本方法的基础上，但是我们这里就不多讲了。可您若发现这是必要的，应该具有相应的背景才可进一步学习。 ---------------------------- 多方面的例子 ------------ 在这个框架中汇集了一些多方面的例子供你们试解答一下。答案提供在框架59，类似于这样的工作有益。 ************************ 偏微分方程 Pg.lp367 框架 58 （续） ------- 1.一张薄型方盘其各边保持温度在-10摄氏度，50摄氏度，100摄氏度和75摄氏度如图所示。。。。？（示意图）？。。。假设稳态温度分布？u（x，y）满足拉普拉斯方程式。。。？（偏微分方程）？。。。，显示这可以由差分方程式逼近。。。？方程式？。。。跨过盘子的任何均匀的网格。所以确定近似值给？u，在已知网格的四个内部节点上。 2. 一个方形盘由线条围住，？x=+-2；y=+-2。盘的温度？西塔o/遵循拉普拉斯方程式。。。？（方程式）？边界温度是：在线。。。？？。。。在线。。。？？。。。在线。。。? ?。。。. 在九个点上寻找温度，这里？i=-1,0,1，j=-1,0,1。 (左手边单元) 3.如果？f（x，y）满足。。。？方程式？。。。，当？f在AB，CD和EF上的值是100时，且是零在PQ和RS上，寻找?f在？O点的近似值。。。。？示意图？。。。（左图所示PQ和RS延伸至无穷双向；CA和DF在左边延伸至无穷。） ****************************************** 偏微分方程 Pg. lp368 框架58（续） ------ 4.如果拉普拉斯方程式由？f（x，y）代替，结果为。。。？（二齐次二维偏微分方程），且已知这为泊松方程式。对于这个方程式对应于？（7.2）公式将是啥？ 5.泊松方程式。。。？二齐次二维偏微分方程式？。。。发生在一定的扭曲问题中。解答这个方程式？0/在网格点如图空心棒所示。它划分为单位方形并且所有的点的值？0/在内和外边界的点都是零。 6. 在OABCDE区域内，?0/满足方程式。。。？（二齐次二维偏微分方程）.？。。。确定到最接近的整数，在网格点，如单元网格所示？0/的值，已知所示在边界的网格点的？0/值。（L.U.）。。。？示意图？。。。 7.寻找解答。。。？（二齐次二维偏微分方程）？。。。，受到条件：？y（0，t）=0，？（1.0，t）=0，？dy（x，0）/dt(偏微分式)=0），？y（x,0）由下表给出 x\|0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 y\|0 0.01 0.02 0.03 0.04 0 -0.04 -0.03 -0.02 -0.01 0 ****************************************** 偏微分方程 Pg.lp369 对于解答，选择一个合适的？k值并进行的足够远以得到运动的整个图像。 8.寻找解答。。。？（偏微分方程）？。。。受到条件：？y（0，t）=0，？y（2，t）=0，？y（x，0）及？偏微分如下表给出 x\|############################### y\|############################### ？偏导y/t\|############################### 连续解答足够远沿时间轴对于一次震荡出现。 9.寻找一根杆子的导热方程的解答。。。？（偏微分方程）？。。。以c2=1/96，取长度1并把它划分为八个间隔。在时间t=0，？西塔0/=100摄氏度通过而不是两端，突然减少到？0摄氏度。假设单元兼容，进行解答直到t=3。 10. 一根米杆延迟沿它的热流满足方程式。。。？（微分方程）？。。。，？c为？1/960.5。在时间t=0，？西塔0/=20度在整个杆子。在x=0流体在400摄氏度被搞得流过端子，而液体在0摄氏度流过另一端。在杆子中间的温度超过35度时时间是多少？H？和K？的值对应的是3000W/m2度及200W/m度，且取h为0.125米。 -------------------------------------------- ****************************************** 单元3-实践数例 Pg.lp372 1. 在一位跳伞运动员从一架飞机跳下之后，紧接着的自由落体运动期间，他的速度，？vm/s在时间？ts已知为？dv/dt=9.81-0.0068v3/2,当t=0时，v=0。运算到小数点5位，使用龙格-库塔方法寻找他在0.4秒之后的速度，采用两步子长度。 2. 使用泰勒序列寻找y，当x=0.1,0.2,0.3 对于方程式。。。？dy/dx=1/6x3+y2,已知y=0.5，当x=0时。那么当x=0.4时，得y，通过米尔尼方法，及当x=0.5时，通过米尔尼修正方法。给出您的结果至4位小数。 3. 当解答线束问题，弧度半径的公式，例如，。。。？{1+（dy/dx）2}3/2/d2y/dx2, 常常取近似为。。。？1/d2y/dx2,假设？（dy/dx）2比1小。那么取微分方程为？EIdy2/dx2=M. 若这种近似不能做，那么微分方程给线束为。。。？EIdy2/dx2=M{1+(dy/dx)2}*3/2 在一个特别的情形？M/EI=10x*2, 所以。。。？（二次微分算式）？。寻找线束的偏斜度，当x=0.2和0.4时，取0.2为步长。假设当x=0时，偏斜度均为零。 4. 一个矩形的盘20厘米（cm）25厘米划分为边长5厘米的小方形。一边保持在100摄氏度，而其它三边在20摄氏度。做到最接近的整数，寻找在网格角的稳态温度。 5.做到最近似的整数，寻找所示内部点的？西塔0/，给出值如图所示。？西塔满足方程式。。。？二次两维度齐次偏微分方程？。。。注意：这里使用的方法与一个举行盘相同。外边不必都与轴平行，可条件是它们通过一系列的节点在已知的？西塔上。。。。？（网格节点示意图）？。。。 **************************************** 单元3-实践例子 Pg.lp373 6. 重做问题5，若现在？西塔满足方程式。。。？（二齐次两维度偏微分方程），小方形的各边是1/2单位。 7. 一根杆子， 0.5米长，被加热，一端在零摄氏度，而另一端在1000摄氏度。经过足够的时间允许整个杆子的温度变得稳定。然后冷端绝缘同时冷却热端，通过允许初始20摄氏度的液体流经它。杆子的温度满足方程式。。。？一二次偏微分方程式？。。。，从冷端测量x，从液体流动开始测量t。取，如文本所述，kc2/h2=1/6，H=3000W/m2摄氏度， ?K=200W/m摄氏度，如果取？h为0.05米，沿着杆子，在t=0.2秒（s）寻找温度分布。全部做到最近的整数。 8. 一根拉伸的绳子在两端固定，（0,0）和（1,0）。然后把绳子的形状做的？y=x/25对于0<=x<=0.5, y=(1-x)/25对于0.5<=x<=1.当在这个位置时，在绳子上的每个点给出一个向着它的平衡位置的速度幅度等于它的位移。从此瞬间测量时间，并取？h=0.1 ，接着计算y直到完成半个周期。对于该特别的运动偏微分方程式是。。。？一二次偏微分方程式？。。。. 9. 给出起始值 x\|0 0.2 0.4 0.6 y\|0 0.0002 0.0032 0.0162 对于方程式。。。？（微分方程）？。。。，使用汉明方法寻找？y，当x=0.8时，且它的修正法寻找？y，当x=1.0。 10. 通过使用射击法，以？h=0.5，寻找至2位小数，y微分的值，当x=0时，即必须取方程式。。。？（二次微分方程）？。。。通过点（0,0）和（1,0）。 11. 使用同时代数方程式法，以h=0.25，寻找到3位小数，？y的值，当?x=0.25,0.5,0.75,对于方程式和条件如问题10。 12. 通过一个泰勒系列的方法，寻找？y的值当x=0.5对于方程式。。。？（二次微分方程）？。。。，已知？y=0和？y’，当x=1.0 及1.5。 13. 已知表？西塔0/（度） \| 70 69 68 67 ？欧米茄（弧度/秒） \|2.083 2.181 2.275 2.365 ？(gsin西塔)/欧米茄l \|5.810 5.513 5.249 5.013 ******************************************** 单元3 --- 实践例子 Pg.lp374-1 以及？g/l=12.88s-2, 解答垂摆方程。。。？（常微分方程）？对于？欧米茄w在？西塔o/=65摄氏度修正到三位小数位置，使用一种高次有限差分方式。提示：亚当-巴西富士公式积分方程式。。。？（方程式）？在区间？A西塔（r，r+1）？是。。。？（算式1）？。。。。。。？（算式2）？。。。（C.E.I.） ******************************* 内容附图页码一致如下：原文附图待续。粤港澳大湾区 2020-8-29 0
2020-8-29 14:31:01　　评论淘帖0 举报相关推荐 • 偏微分方程式数值解之计算机Fortran算法详解附图（续ii） 1733 • 偏微分方程式数值解之计算机Fortran算法详解附图（续i） 2029 • 算法大全_偏微分方程的数值解 0 • 微分方程式的建立与求解 0 • 量子技术使用线性方程式作为解非线性微分方程式 1628 • 计算机Fortran编程实现偏微分方程拉普拉斯变换 2976 • 异性四阶偏微分方程耦合二阶偏微分方程的图像放大算法 0 • 差分方程与原偏微分方程的区别与特点 1648 • 有限差分法求解偏微分方程 1930 • 基于偏微分方程的闪光照相图像修补算法 22