fortran 90 关于线性病态方程组解答程序
fortran数值计算基础
数值计算基础 目录 实验一 直接法解线性方程组的 .......................................................... 2 实验二 插值方法 ............................................................................. 12 实验三 数值积分 ............................................................................... 6 实验四 常微分方程的数值解 .............................................................. 8 实验五 迭代法解线性方程组与非线性方程 ...................................... 10 实验一 直接法解线性方程组 一、实验目的 掌握全选主元消去法与高斯-塞德尔法解线性方程组。 二、实验内容 分别写出Guass列选主元消去法与追赶法的算法，编写程序上机调试出结果，要求所编程序适用于任何一解线性方程组问题，即能解决这一类问题，而不是某一个问题。实验中以下列数据验证程序的正确性。 1、用Guass列选主元消去法求解方程组
5.58.37.33.47.11.85.16.93.51.53.25.2321xxx
2、用追赶法求解方程组
0000102100002100002100002100002
54321
xxxx
x
三、实验仪器设备与材料 主流微型计算机 四、实验原理 1、Guass列选主元消去法 对于AX =B 1）、消元过程：将（A|B）进行变换为)~|~(BA，其中A~是上三角矩阵。即：
数值分析（Hilbert矩阵）病态线性方程组的求解Matlab程序
数值分析（Hilbert矩阵）病态线性方程组的求解Matlab程序（Hilbert 矩阵）病态线性方程组的求解理论分析表明，数值求解病态线性方程组很困难。
考虑求解如下的线性方程组的求解Hx = b ，期中H 是Hilbert 矩阵，()ij n n H h ?=，11ij h i j =+-，i ，j = 1,2，…，n 1. 估计矩阵的2条件数和阶数的关系 2. 对不同的n ，取(1,1,,1)nx =∈K ?，分别用Gauss 消去，Jacobi 迭代，Gauss-seidel 迭代，SOR 迭代和共轭梯度法求解，比较结果。
3. 结合计算结果，试讨论病态线性方程组的求解。
第1小题：condition.m %第1小题程序t1=20;%阶数n=20x1=1:t1;y1=1:t1;for i=1:t1H=hilb(i);y1(i)=log(cond(H));endplot(x1,y1);xlabel('阶数n');ylabel('2-条件数的对数（log(cond(H)）');title('2-条件数的对数（log(cond(H)）与阶数n 的关系图');t2=200;%阶数n=200x2=1:t2;y2=1:t2;for i=1:t2H=hilb(i);y2(i)=log(cond(H));endplot(x2,y2);xlabel('阶数n');ylabel('2-条件数的对数（log(cond(H)）');title('2-条件数的对数（log(cond(H)）与阶数n 的关系图');画出Hilbert 矩阵2-条件数的对数和阶数的关系n=200时n=20时从图中可以看出，1）在n小于等于13之前，图像近似直线log（cond（H））~1.519n-1.8332）在n大于13之后，图像趋于平缓，并在一定范围内上下波动，同时随着n的增加稍有上升的趋势第2小题：solve.m%m第2小题主程序N=4000;xGauss=zeros(N,1);xJacobi=zeros(N,1);xnJ=zeros(N,1);xGS=zeros(N,1);xnGS=zeros(N,1);xSOR=zeros(N,1);xnSOR=zeros(N,1);xCG=zeros(N,1);xnCG=zeros(N,1);for n=1:N;x=ones(n,1);t=1.1;%初始值偏差x0=t*x;%迭代初始值e=1.0e-8;%给定的误差A=hilb(n);b=A*x;max=100000000000;%可能最大的迭代次数w=0.5;%SOR迭代的松弛因子G=Gauss(A,b);[J,nJ]=Jacobi(A,b,x0,e,max); [GS,nGS]=G_S(A,b,x0,e,max);[S_R,nS_R]=SOR(A,b,x0,e,max,w);[C_G,nC_G]=CG(A,b,x0,e,max);normG=norm(G'-x);xGauss(n)=normG;normJ=norm(J-x);nJ;xJacobi(n)=normJ;xnJ(n)=nJ;normGS=norm(GS-x);nGS;xGS(n)=normGS;xnGS(n)=nGS;normS_R=norm(S_R-x); nS_R;xSOR(n)=normS_R; xnSOR(n)=nS_R; normC_G=norm(C_G-x); nC_G;xCG(n)=normC_G; xnCG(n)=nC_G;endGauss.m%Gauss消去法function x=Gauss(A,b)n=length(b);l=zeros(n,n);x=zeros(1,n); %消去过程for i=1:n-1for j=i+1:nl(j,i)=A(j,i)/A(i,i);for k=i:nA(j,k)=A(j,k)-l(j,i)*A(i,k); endb(j)=b(j)-l(j,i)*b(i);endend%回代过程x(n)=b(n)/A(n,n);for i=n-1:-1:1c=A(i,:).*x;x(i)=(b(i)-sum(c(i+1:n)))/A(i,i);endJacobi.m%Jacobi迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m 可能最大的迭代次数function [x,n]=Jacobi(A,b,x0,e,m)n=length(A);D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B=D\(L+U);f=D\b;x=B*x0+f;n=1;while norm(x-x0)>ex0=x;x=B*x0+f;n=n+1;if n>mdisp('Jacobi迭代次数过多，迭代可能不收敛');break;endendG_S.m%Gauss-Seidel迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数function [x,n]=G_S(A,b,x0,e,m)n=length(A);D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B=(D-L)\U;f=(D-L)\b;x=B*x0+f;n=1;while norm(x-x0)>ex0=x;x=B*x0+f;n=n+1;if n>mdisp('Gauss-Seidel迭代次数过多，迭代可能不收敛');break;endendSOR.m%SOR超松弛迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数,w松弛因子function [x,n]=SOR(A,b,x0,e,m,w)n=length(A);D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B=(D-w*L)\((1-w)*D+w*U);f=(D-w*L)\b*w;x=B*x0+f;n=1;while norm(x-x0)>ex0=x;x=B*x0+f;n=n+1;if n>mdisp('SOR超松弛迭代次数过多，迭代可能不收敛');break;endendCG.m%CG共轭梯度法,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数function [x,n]=CG(A,b,x0,e,m)r=b-A*x0;p=r;alpha=(r'*r)/(p'*(A*p));x=x0+alpha*p;r1=b-A*x;n=1;while norm(r1)>ebelta=(r1'*r1)/(r'*r);p=r1+belta*p;r=r1;x0=x;alpha=(r'*r)/(p'*(A*p));x=x0+alpha*p;r1=b-A*x;n=n+1;if n>mdisp('CG共轭梯度法迭代次数过多，迭代可能不收敛');break;endend。
FORTRAN90指导
FORTRAN90指导辅导资料――*****90 楚红*****90学习指导一、*****90基础知识1．源程序及其构成书写格式：主程序定义语句、结束语句注意：在没有程序名称的时候，程序开头的*****也不要出现。
如果使用了程序名称且在END语句中出现，则END语句中的*****不能省略。
语句行可以是0~132个字符；除赋值语句外，每个语句都要使用关键字开头。
如果希望一行中出现多个语句，语句一定要用分号隔开。
空格不能随便使用，关键字、变量和常量名以及操作符中的空格会使字符失去其原有的含义。
但它们之间一定要加空格。
注释行：以感叹号为标记，或“C”、“*”（*****90中不提倡此用法）续行标记：在句末尾添加续行符，如果将关键字分成两行，则下一行开头也要加续行标记。
如果续行符出现在注释语句中，则失去了续行的功能。
2．*****90字符集1 26个英文字母（大小写字母等价）2 10个阿拉伯数字3 下划线4 21个特殊字符：空格= + - * / （），. ‘ : ！“ % ；? $5 其他字符（只可以出现在字符常量、字符串编辑描述符、注释和输入输出记录中）3．基本数据类型1 整型KIND值可以为1、2、4 类型说明关键字：*****KIND值决定数据的范围。
对于整数没有误差。
I=3/2=1；I=1/2=0 *****（KIND=4）：：A ******4：：A *****（4）：：A2 实型KIND值可以为4（单精度7位有效数字）、8（双精度15~17位有效数字）类型说明关键字：REAL小数表示形式、指数表示形式（规格化的指数形式）错误表示形式：E34、.E34、0.14E2.3 8开3次方：8.0**（1.0/3.0）3 复型KIND值为4、8是实数的有序对，将两个实数中间用逗号分隔，然后放在一对括号中类型说明关键字：***** 复数的运算：+、―、*、/4 字符型KIND值为国家语言种类由一对单撇号或一对双撇号之间的字符序列组成。
病态希尔伯特方程组解法-高等数值分析大作业
病态方程组的求解理论分析表明，数值求解病态线性方程组很困难，考虑求解如下的线性方程组的求解：Hx=b ，其中H 是Hilbert 矩阵，H=(hij)n×n ，h ij =1i+j−1,ii ,jj =1,2…..nn1.估计矩阵的2-条件数和矩阵阶数的关系；2.对不同的n ，取x =(1,1,1…,1)∈R n ，利用Hx=b 得出b.然后分别用Guass 消去，Jacobi 迭代，Gauss-seidel 迭代，SOR 迭代求解方程： Ax=b ；把近似解x (k )与准确解x 对比，计算结果；3.综合计算结果，讨论病态方程组的求解。
解：1.用MATLAB 计算1-9阶Hilbert 矩阵的2-条件数得到的结果如下表所示（cond2（Hnn ）表示n 阶Hilbert 矩阵的2-条件数）：对表中数据进行分析可发现，随矩阵的阶数的增加，2-条件数越来越大，后项与前项之比约为30，2-条件数呈近似指数式增长。
对cond2（H nn ）取以10为底的对数，可发现在前几项log10(cond2(Hn))随n 的增加约呈线性增加关系。
n 为20时，log10(cond2(Hn))与n 的关系与下图所示：此时可发现当n ≤13时，log10(cond2(Hn))与n 接近于线性关系，n >13后，log10(cond2(Hn))变化十分缓慢。
为了对图像的线性进行比较，对n 为1-13部分1-9阶hilbert 矩阵的2-条件数n 12 3 4 5 6 7 8 9 cond2(H nn ) 119.2815 524.0568 1.55E+04 4.77E+05 1.49E+07 4.75E+08 1.53E+10 4.93E+11进行线性拟合，结果如下图所示：从上图可看出，n较小时，log10(cond2(Hn))与n成近似线性关系。
继续增大n，n=200时，log10(cond2(Hn))与n如下图所示：由图可看出，n为200时，log10(cond2(Hn))与n的关系与n=20时相近，在n较小时log10(cond2(Hn))与n的关系接近线性关系，当n大于一定的值后log10(cond2(Hn))变化十分缓慢，但有上升的趋势。
fortran课后习题答案
fortran课后习题答案Fortran课后习题答案Fortran是一种编程语言，广泛应用于科学计算和工程领域。
学习Fortran的过程中，课后习题是非常重要的一部分，通过解答这些习题，可以巩固所学的知识，并且提升自己的编程能力。
在本文中，我将为大家提供一些Fortran课后习题的答案，希望对大家的学习有所帮助。
1. 编写一个Fortran程序，计算并输出1到100之间所有偶数的和。
```fortranprogram sum_even_numbersimplicit noneinteger :: i, sumsum = 0do i = 1, 100if (mod(i, 2) == 0) thensum = sum + iend ifend doprint *, "The sum of even numbers from 1 to 100 is", sumend program sum_even_numbers```2. 编写一个Fortran程序，计算并输出用户输入的两个数的乘积。
```fortranprogram multiply_numbersimplicit nonereal :: num1, num2, productprint *, "Enter the first number:"read *, num1print *, "Enter the second number:"read *, num2product = num1 * num2print *, "The product of", num1, "and", num2, "is", productend program multiply_numbers```3. 编写一个Fortran程序，计算并输出用户输入的一个数的平方根。
求解非线性方程组的几种方法及程序实现
求解非线性方程组的几种方法及程序实现
求解非线性方程组一直是理论数学和应用数学研究的重点，并采用不同的方法得到准确的结果。
它们可以分为几种类型：
1. 用以绘图的方法解非线性方程组：该方法充分利用结合几何和数理的原理，给出非线性方程组的解，而不用对系数的解的表达式求解手段。
主要是利用可绘图的几何空间分析，它可以帮助理解问题本身，还可以很容易看出非线性方程组的解。
2. 用迭代法求解非线性方程组：这是一种常用的方法，它通过不断迭代收敛求解非线性方程组。
基本思想是通过构造一个迭代函数，其初始值和原始非线性方程组尽可能接近，然后不断迭代收敛求解非线性方程组。
3. 用强调法求解非线性方程系统：这是基于梯度的一种方法，它利用一个概念，即局部线性化，可以降低维数、转化为一个拐点，最后强化搜索全局解。
4. 用牛顿-拉夫逊方法求解非线性方程组：这是一种准确、快速的非线性方程组求解方法，主要利用牛顿迭代法搜索解的收敛性，加上一些拉夫逊的加速策略得到最终的结果。
5. 用幂法求解非线性方程组：幂法也称为指数序列，是一种重要的求解非线性方程组的方法，基本原理是利用指数的累加和误差的减少，从而最终得到非线性方程组的解。
6. 用逐步逼近法求解非线性方程组：逐步逼近法也称为分步变程法，是一种用于求解非线性方程组的简单方法，其基本思想是用不同的参数，在给定的范围内，逐步逼近目标解。
这些方法的程序实现略有不同，可以利用编程语言比如C、Fortran、Python等，编写程序完成求解。
可以采用函数求解、循环求解、行列式求解或者混合的算法等不同的方式实现，甚至可以用深度学习方法求解有些复杂的非线性方程组。
fortran90整理
fortran90整理1语句编译1.Build—Compile：编译；Build—Build：连接；Build—Exetuce：运⾏；或单击⼯具栏相应按钮。
注意：a、保存⽂件时将⾃动创建同名的project⽂件，形成*.dsp⽂件；b、同时还将⾃动创建同名的workspace，形成*.opt和*.dsw⽂件；c、编译连接后⾃动形成Debug⽬录，该⽬录中存放编译连接后⽂件。
如：*.obj,*.lnk,*.exe等2.Free Format（⾃由格式）1.！感叹号后⾯的⽂本都是注释。
2.每⾏可以编写132个字符。
3.⾏号放在每⾏程序的最前⾯。
4.⼀⾏程序代码的最后如果是符号&，代表下⼀⾏程序会和这⼀⾏连接。
如果⼀⾏程序代码的开头是符号&，代表它会和上⼀⾏程序连接。
3.书写格式⾏的书写(⾏的长度、分⾏、续⾏)1⼀⾏可以是0~132个字符，空格有意义，2语句最长不超过2640个字符3⼀⾏可以有多个语句，⽤“;”分隔4⼀个语句可分⾏写，读⾏标记为&（放在尾部），但如为关键字，⾸尾均加&。
5最多可有511个续⾏。
4.语句的分类注释语句：！后的所有字符都被编译器忽略1.可独占⼀⾏，可在其它语句之后，a)空⾏为注释⾏（固定格式⽤C和*）2.说明语句：⽤于说明变量的类型、属性等3.可执⾏语句：输⼊、赋值、输出……5.语句有位置规定：说明语句必须出现在可执⾏语句之前，格式说明语句（FORMAT语句）除外。
6.标志符⼩结注释标志符：1⾃由格式：！固定格式：C *2语句分隔符：分号；（仅⾃由格式可以使⽤）3续⾏符：⾃由格式：&4申明标号：1到5位⽆符号整数5空格：关键字、变量、常量内部不能⽤空格，但相邻两者之间须⽤空格6.FORTRA90源程序基本结构1、FORTRAN90程序是⼀种分块结构，由若⼲个程序单元块组成：主程序、外部⼦程序、模块、块数据单元⽆论是主程序单元，还是⼦程序单元，都是独⽴的程序单位，应该独⽴编写，它们的形式相似。
矩阵病态线性代数方程组的求解
实验一病态线性代数方程组的求解1.估计Hilbert矩阵2-条件数与阶数的关系运行tiaojianshu.m输入m=10 可以得到如下表的结果2.选择不同维数，分别用Guass消去（LU分解），Jacobi迭代，GS 迭代，SOR迭代求解，比较结果。
说明：Hx=b，H矩阵可以由matlab直接给出，为了设定参考解，我们先设x为分量全1的向量，求出b,然后将H和b作为已知量，求x，与设定的参考解对比。
对于Jacobi迭代，GS迭代，SOR迭代，取迭代初值x0为0向量，迭代精度eps=1.0e-6，迭代次数<100000, SOR迭代中w=1.2和0.8分别计算。
取不同的n值，得到如下结果:对于Guass法，可以看出来，随着n的增大，求解结果误差变大，这是因为随着n增大，系数矩阵的条件数变大，微小的扰动就容易造成很大的误差。
最后得不到精确解。
对于Jacobi迭代，计算结果为Inf，说明是发散的。
对于GS迭代和SOR迭代，结果是收敛的，但是可以看出迭代次数比较多，并且对于不同维数GS和SOR收敛速度不一样，有时候GS快，有时SOR快。
对SOR取不同的w迭代速度也不一样，存在一个最优的松弛因子w。
并且可以知道，迭代次数多少跟初值x0也有关系。
3.讨论病态问题求解的算法。
通过上面的实验分析，可以看出，求解病态矩阵的时候要小心，否则可能得不到所要求的精确度。
可以采用高精度运算，用双倍多倍字长，使得由于误差放大而损失若干有效数字位之后，还能保留一些有效位。
另外可以通过对原方程作某些预处理，降低系数矩阵的条件数，因为cond(aA)=cond(A),所以不能通过将每一个方程乘上相同的常数来达到这个目标，可考虑将矩阵的每一行和每一列分别乘上不同的常数，亦即找到可逆的对角阵D1和D2将方程组化为D1AD2y=D1b，x=D2y这称为矩阵的平衡问题，但是这样计算量比原问题本身要多。
或者通过变分原理将求解线性方程组的问题转化为等价的求解无约束函数最优化问题的极小值等等，可以参考[1]郑洲顺，黄光辉，杨晓辉求解病态线性方程组的混合算法实验一所编程序如下：1.求条件数tiaojianshu.mm=input ('input m:=') ;N=[1:m];for i=1:mn=N(i);H=hilb(n);k=cond(H);disp('矩阵的阶数')disp(n)disp('矩阵')disp(H)disp('矩阵的条件数')disp(k)end2.Guass法①Guass.mfunction guass(n)n=input('请输入系数矩阵的维数n=');H=hilb(n);x_exact=ones(n,1);b=H*x_exact;x=Doolittle(H,b)a=input('是否继续,继续请按1,结束请按0:');if a==1guass(n);else end②Doolittle.mfunction x=Doolittle(A,b)% LU分解求Ax=b 的解N=size(A);n=N(1);L=eye(n,n);%L的对角元素为1U=zeros(n,n); %U为零矩阵U(1,1:n)=A(1,1:n)%U第一行L(1:n,1)=A(1:n,1)/U(1,1)%L第一列for k=2:nfor i=k:nU(k,i)=A(k,i)-L(k,1:(k-1))*U(1:(k-1),i)endfor j=(k+1):nL(j,k)=(A(j,k)-L(j,1:(k-1))*U(1:(k-1),k))/U(k,k) endendy=solveDownTriangle(L,b);%调用下三角系数矩阵求解函数x=solveUpTriangle(U,y);%调用上三角系数矩阵求解函数③solveDownTriangle.mfunction x=solveDownTriangle(A,b)%求下三角系数矩阵的线性方程组Ax=bN=size(A);n=N(1);for i=1:nif (i>1)s=A(i,1:(i-1))*x(1:(i-1),1);elses=0;endx(i,1)=(b(i)-s)/A(i,i);end④solveUpTriangle.mfunction x=solveUpTriangle(A,b)% 求上三角系数矩阵的线性方程组Ax=bN=size(A);n=N(1);for i=n:-1:1if (i<n)s=A(i,(i+1):n)*x((i+1):n,1);elses=0;endx(i,1)=(b(i)-s)/A(i,i);end3.Jacobi法function [x,n]=jacobi(a,x0)a=input('请输入系数矩阵的维数n=');x0=input('请输入迭代初始向量x0');eps=1.0e-6;%解的精度控制m=10000;%迭代步数控制H=hilb(a);%生成h矩阵x_exact=ones(a,1);%求出x精确值b=H*x_exact;D=diag(diag(H));L=-tril(H,-1);U=-triu(H,1);B=D\(L+U);f=D\b;x=x0;n=0;tol=1;while tol>=epsx=B*x0+f;n=n+1;tol=norm(x-x0);x0=x;if (n>=m)disp('迭代次数太多，可能不收敛');return;endenddisp(x)disp(n)4.GS法function [x,n]=gauseidel(a,x0)a=input('请输入系数矩阵的维数n=');x0=input('请输入迭代初始向量x0');eps=1.0e-6;%解的精度控制m=100000;%迭代步数控制H=hilb(a);%生成h矩阵x_exact=ones(a,1);%求出x精确值b=H*x_exact;D=diag(diag(H));L=-tril(H,-1);U=-triu(H,1);G=(D-L)\U;f=(D-L)\b;x=x0;n=0;tol=1;while tol>=epsx=G*x0+f;n=n+1;tol=norm(x-x0);x0=x;if (n>=m)disp('迭代次数太多，可能不收敛');return;endenddisp(x)disp(n)5.SOR法function [x,n]=SOR(a,x0,w)a=input('请输入系数矩阵的维数n=');x0=input('请输入迭代初始向量x0=');w=input('请输入w=');eps=1.0e-6;%解的精度控制m=100000;%迭代步数控制H=hilb(a);%生成h矩阵x_exact=ones(a,1);%求出x精确值b=H*x_exact;if (w<=0||w>=2)error;return;endD=diag(diag(H));L=-tril(H,-1);U=-triu(H,1);B=inv(D-L*w)*((1-w)*D+w*U);f=w*inv((D-L*w))*b;x=x0;n=0;tol=1;while tol>=epsx=B*x0+f;n=n+1;tol=norm(x-x0);x0=x;if (n>=m)disp('迭代次数太多，可能不收敛');return;endenddisp(x)disp(n)。
FORTRAN_90标准函数
附录 FORTRAN 90标准函数符号约定：
●I代表整型;R代表实型;C代表复型;CH代表字符型;S代表字符串;L代表逻辑型;A代表数组;P代表指针;T代表派生类型;AT为任意类型。
●s:P表示s类型为P类型(任意kind值)。
s:P(k)表示s类型为P类型(kind值=k)。
●[…]表示可选参数。
●*表示常用函数。
表1 数值和类型转换函数
FORTRAN 90程序设计
表2 三角函数
注：三角函数名前有C、D的函数为复数、双精度型函数。
表3 指数、平方根和对数函数
-394-
附录 FORTRAN 90标准函数简表
表4 参数查询函数
表5 实数检测和控制函数
表6 字符处理函数
395
FORTRAN 90程序设计
表7 二进制位操作函数
表8 数组运算、查询和处理函数
-396-
附录 FORTRAN 90标准函数简表
注: 参数m指逻辑型掩码数组,指明允许操作的数组元素。
缺省掩码数组指对数组所有元素进行操作。
397。
(最新整理)FORTRAN90第五章过程
2021/7/26
26
例:编写程序计算一个角度的正切和余切值. Function triangle(x,f1,f2)result(tria_result) real x, tria_result tria_result=f1(x)/f2(x) End Function triangle
Program main
End function F
End 2021/7/26
15
§5.3 子例行程序
§5.3.1 外部子例行程序
子例行程序的定义: Subrotine 子例行程序名([虚参表]) …过程体… ＥＮＤ [Subrotine [子例行程序名]]
如果没有参数,则定义格式中()可有可无。
子例行程序调用方法:
ＣＡＬＬ 子例行程序名(实参表)
i1 k
在x=5.0,5.2,…k61.8处,计算对应的一次至五
次多项式
function Fn(n,x) real x;integer n fn = 0;Xi=1.0 Do I = 1,n S = sum(I) Xi = Xi*X Fn = fn+XI/s enddo end 2021/7/26
Function sum(k) Sum=0.0 Do j=1,k Sum=sum+j End do end
宿主程序的说明语句与数据定义,是全局数 据,内部过程中都可以使用,并且内部过程对 这些数据的赋值对宿主程序也有效。
当局部数据与全局数据同名时,局部数据屏 蔽全局数据。在内部过程中对变量的使用与 赋值都只与局部变量有关。
内部过程必须放在宿主程序的最后。宿主
程序只有一条END语句在内部过程的后
面。 2021/7/26
…函数体… END [Function [fun]]
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
program lichunxin parameter (n=3,np=3) dimension a(n,n),indx(n),b1(n),b2(n),x1(n),x2(n),aa(n,n),c(n),al(n,n) real cond integer n real z a(1,1)=1.0/2;a(1,2)=1.0/3;a(1,3)=1.0/4 a(2,1)=1.0/3;a(2,2)=1.0/4;a(2,3)=1.0/5 a(3,1)=1.0/4;a(3,2)=1.0/5;a(3,3)=1.0/6 b1=(/.95,.67,.52/) do i=1,n m=1 do i1=i+1,i+n a(i,m)=1.0/i1 m=m+1 end do end do open(10,file='f:\a1.txt') write(10,100)a write(10,100)b1 100 format(3f15.6) close(10) open(10,file='f:\a1.txt') read(10,100)a read(10,100)b1 close(10) do i=1,n x1(i)=b1(i) do j=1,n aa(i,j)=a(i,j) enddo enddo call LUD(aa,n,np,indx,d) call LUB(aa,n,np,indx,x1) write(*,'(/1x,A)')"输出当b3=0.52时方程组的解：" write(*,10)(x1(i),i=1,n) 10 format(1x,3f15.6) call MPROVE(a,aa,n,np,indx,b1,x1) write(*,'(/1x,A)')"输出当b3=0.52时方程组改善后的解：" write(*,30)(x1(i),i=1,n) 30 format(1x,3f15.6) do l=1,n c(l)=0. do j=1,n c(l)=c(l)+a(l,j)*x1(j) enddo enddo print 400,"输出原矩阵验证结果：" 400 format(/1x,A) write(*,'(1x,3f15.6//)')(c(l),l=1,n) a(1,1)=1./2;a(1,2)=1./3;a(1,3)=1./4 a(2,1)=1./3;a(2,2)=1./4;a(2,3)=1./5 a(3,1)=1./4;a(3,2)=1./5;a(3,3)=1./6 b2=(/.95,.67,.53/) do i=1,n m=1 do i1=i+1,i+n a(i,m)=1.0/i1 m=m+1 end do end do open(11,file='f:\a2.txt') write(11,100)a write(11,101)b2 101 format(3f15.6) close(11) open(11,file='f:\a2.txt') read(11,100)a read(11,100)b2 close(11) do i=1,n x2(i)=b2(i) do j=1,n aa(i,j)=a(i,j) enddo enddo call LUD(aa,n,np,indx,d) call LUB(aa,n,np,indx,x2) write(*,'(/1x,A)')"输出当b3=0.53时方程组解：" write(*,40)(x2(i),i=1,n) 40 format(1x,3f15.6) call MPROVE(a,aa,n,np,indx,b2,x2) write(*,'(/1x,A)')"输出当b3=0.53时方程组改善后的解：" write(*,60)(x2(i),i=1,n) 60 format(1x,3f15.6) do l=1,n c(l)=0. do j=1,n c(l)=c(l)+a(l,j)*x2(j) enddo enddo print 200,"输出更改矩阵验证结果：" 200 format(/1x,A) write(*,'(1x,3f15.6)')(c(l),l=1,n) open(11,file='f:\a2.txt') read(11,100)a close(11) call fanshu(a,n,fs1) print*,"矩阵A的范数:""fs1=",fs1 call sjz(a,n) sum=1.0 do i=1,n sum=sum*a(i,i) end do z=sum if(z==0)then print*,'系数矩阵为退化矩阵' else if(z/=0)then print*,'此时非退化行列式值为' print*,z end if open(11,file='f:\a2.txt') read(11,100)a read(11,100)b2 close(11) call njz(a,n,z,al) call fanshu(al,n,fs3) fs2=fs3 print*,"A的逆矩阵的范数:""fs2=",fs2 print 300,"条件数cond=fs1*fs2=" 300 format(/1x,A) cond=fs1*fs2 write(*,70) cond 70 format(1x,f15.6) end
subroutine lud(a,n,np,indx,d) parameter (nmax=100,tiny=1.0e-20) dimension a(np,np),indx(n),v(nmax) real*8 aamax d=1.0 do i=1,n aamax=0.0 do j=1,n if(abs(a(i,j))>aamax) aamax=abs(a(i,j)) enddo if (aamax==0.) pause 'singular matrix' v(i)=1.0/aamax enddo do j=1,n if (j>1) then do i=1,j-1 sum=a(i,j) if(i>1) then do k=1,i-1 sum=sum-a(i,k)*a(k,j) enddo a(i,j)=sum endif enddo endif aamax=0.0 do i=j,n sum=a(i,j) if (j>1) then do k=1,j-1 sum=sum-a(i,k)*a(k,j) enddo a(i,j)=sum endif dum=v(i)*abs(sum) if (dum>=aamax) then imax=i aamax=dum endif enddo if(j/=imax) then do k=1,n dum=a(imax,k) a(imax,k)=a(j,k) a(j,k)=dum enddo d=-d v(imax)=v(j) endif indx(j)=imax if(j/=n) then if (a(j,j)==0.0) a(j,j)=tiny dum=1.0/a(j,j) do i=j+1,n a(i,j)=a(i,j)*dum enddo endif enddo if(a(n,n)==0) a(n,n)=tiny end subroutine LUD
subroutine LUB(a,n,np,indx,b) dimension a(np,np),indx(n),b(n) ii=0 do i=1,n ll=indx(i) sum=b(ll) b(ll)=b(i) if(ii/=0) then do j=ii,i-1 sum=sum-a(i,j)*b(j) enddo else if (sum/=0.0)then ii=i endif b(i)=sum enddo do i=n,1,-1 sum=b(i) if(ido j=i+1,n sum=sum-a(i,j)*b(j) enddo endif b(i)=sum/a(i,i) enddo end subroutine LUB
subroutine MPROVE(a,alud,n,np,indx,b,x) parameter (nmax=100) real a(np,np),alud(np,np),b(n),x(n),r(nmax) integer i,n,indx(np) real*8 sdp do i=1,n sdp=-b(i) do j=1,n sdp=sdp+dble(a(i,j))*dble(x(j)) enddo r(i)=sdp enddo call LUB(alud,n,np,indx,r) do i=1,n x(i)=x(i)-r(i) enddo end subroutine MPROVE
subroutine fanshu(b,n,fs) dimension b(n,n),c(n) do i=1,n c(i)=0.0 do j=1,n c(i)=c(i)+abs(b(i,j)) enddo enddo fs=0.0 do i=1,n if(fsfs=c(i) endif enddo fs=fs end
subroutine njz(a,n,z,c) real a(n,n),g(n-1,n-1),d1((n-1)*(n-1)),c(n,n) integer n real sum,z do i=1,n do j=1,n m=1 do k1=1,n do k=1,n if(k/=i.and.k1/=j)then d1(m)=a(k,k1) m=m+1 end if end do end do