fortran线性病态方程组问题

合集下载

范德蒙矩阵形式下的病态线性方程组求解

由表 1 的数据可知，随着范德蒙德矩阵阶数的增加，其 2- 条件数也越来越大，病态性也越来越严重参见文献[6-7]。为更直观地了解阶数与条件
基金项目：2018 年长治学院课题“非线性算子的不动点定理及其应用”（ZC201811）；长治学院“1331 工程”人才培养质量提升计划项目（200628）
n
移 bi= aij，i=1，2，…n j=1
下面用新主元加权迭代法对这个线性方程组进行求解，得到的结果如表 3 所示。
表 3 解的近似值（n=10 加权因子为籽=1.000001）
近似值 x1 x2 x3 x4 x5 x6 x7 x8 x9 x10
单参数迭代法 0.999999 1 0.999997 1 0.999999 1 1 1 1 1
表 4 加权因子与绝对误差
文章选取以系数矩阵为范德蒙德矩阵的病态线性方程组，借鉴文献[4]提供的单参数迭代法和文献[5]中的新主元加权迭代法，对病态线性方程组进行分析求解。结果表明选取的迭代方法切实可行，对分析此病态线性方程组有很大帮助。
1 范德蒙德矩阵的病态性分析
选取 n 阶的范德蒙德矩阵如下：
1 杉山
山
1
12
…
1 山
数的增加，2- 条件数越来越大，病态程度也越来越严重。然后选用单参数迭代法和新主元加权迭代法分别
对以系数矩阵是范德蒙德矩阵的病态线性方程组进行求解，从数值结果可以看出，选取适当的加权因子对
求解此病态线性方程组同样有较好的精度和收敛速度。
关键词：病态线性方程组；范德蒙矩阵；单参数迭代法；新主元加权迭代法
2019 年 4 月第 36 卷第 2 期
长治学院学报 Journal of Changzhi University

fortran求解方程复数结果

fortran求解方程复数结果Fortran不直接支持复数运算，但可以通过引入虚数单位来实现复数运算。

以下为一个Fortran程序示例，演示了如何使用复数运算库来求解方程并输出复数结果。

```fortranprogram complex_equationimplicit noneinteger, parameter :: dp = selected_real_kind(15)complex(dp) :: x, solution! 引入复数运算库use iso_fortran_env! 定义方程和初始猜测solution = cmplx(1.0_dp, 1.0_dp, kind=dp)! 使用复数运算库进行迭代求解do while (.true.)x = newton_method(solution, dp)if (abs(x - solution) < 1e-6_dp) exitsolution = xend do! 输出结果write(*, *) "The solution is ", xcontainsfunction newton_method(x, dp) result(new_x)implicit nonecomplex(dp), intent(in) :: xcomplex(dp) :: new_x, f, df! 定义方程和其导数f = x**2 - (3.0_dp, 4.0_dp)df = 2.0_dp * x! 使用牛顿法更新猜测new_x = x - f / dfend function newton_methodend program complex_equation```在上述示例中，我们定义了一个名为`newton_method`的函数，它使用牛顿法进行迭代求解方程。

然后，我们使用初始猜测和复数运算库来求解方程，并输出结果。

请注意，该示例仅仅在复数部分存在的情况下展示了Fortran的求解方程示例，如果方程具有实数解，可以在实数范围内解决该问题。

fortran求解方程组

fortran求解方程组在Fortran中，可以使用数值方法来求解方程组。

一种常用的数值方法是迭代法，其中最常见的方法是牛顿-拉夫逊方法（Newton-Raphson method）。

以下是一个使用Fortran实现牛顿-拉夫逊方法求解方程组的示例：```fortranprogram newton_raphsonimplicit noneinteger, parameter :: n = 2 ! 方程组中的未知数个数integer, parameter :: max_iter = 100 ! 最大迭代次数real, parameter :: tolerance = 1e-6 ! 精度要求real :: x(n) ! 未知数向量real :: f(n) ! 方程组函数值向量real :: J(n,n) ! 方程组雅可比矩阵real :: dx(n) ! 迭代步长向量integer :: iter ! 迭代次数! 初始化未知数向量x = [1.0, 1.0]iter = 0 ! 初始化迭代次数do while (iter < max_iter)iter = iter + 1! 计算方程组函数值向量和雅可比矩阵call calculate_f(x, f)call calculate_J(x, J)! 解线性方程组 J * dx = -f，计算迭代步长 call solve_linear_equation(J, f, dx)! 更新未知数向量x = x + dx! 判断迭代是否收敛if (maxval(abs(f)) < tolerance) thenexitend ifend doif (iter >= max_iter) thenprint *, "迭代失败：达到最大迭代次数" elseprint *, "迭代成功：达到收敛精度"end ifprint *, "迭代次数：", iterprint *, "解向量：", xcontains! 计算方程组函数值向量subroutine calculate_f(x, f)real, intent(in) :: x(n)real, intent(out) :: f(n)! 实现方程组函数值的计算，将结果存储到 f 中end subroutine calculate_f! 计算方程组雅可比矩阵subroutine calculate_J(x, J)real, intent(in) :: x(n)real, intent(out) :: J(n,n)! 实现方程组雅可比矩阵的计算，将结果存储到 J 中end subroutine calculate_J! 解线性方程组subroutine solve_linear_equation(A, b, x)real, intent(in) :: A(n,n)real, intent(in) :: b(n)real, intent(out) :: x(n)! 实现线性方程组的求解，将结果存储到 x 中end subroutine solve_linear_equationend program newton_raphson```在上面的示例中，你需要根据具体的方程组函数和雅可比矩阵的计算方法，实现相应的子程序`calculate_f`、`calculate_J` 和`solve_linear_equation`。

计算机程序设计基础—FORTRAN实验设计报告线性方程组求解问题

中南大学本科生课程设计(实践)任务书、设计报告（计算机程序设计基础—FORTRAN）题目线性方程组求解问题学生姓名陈晨指导教师刘胤宏学院土木工程学院专业班级土建类工程试验班学生学号18计算机基础教学实验中心2012年 6 月29日Fortran 课程设计实验报告之线性方程组求解问题题目重现：一物理系统可用下列线性方程组来表示：⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡θ-θ-θθθ--θg m g m N N a a m m m m 2121212111001cos 00sin 00cos 0sin 0sin cos 从文件中读入m1、m2和θ的值，求a1、a2、N1 和N2的值。

其中g 取，输入θ时以角度为单位。

要求如下：（1）分别用两种方法（例如高斯消去法、矩阵求逆法、三角分解法、追赶法等），定义求解线性方程组Ax=b 的子程序，要求该子程序能求解任意线性方程组。

（2）在主程序中分别调用上面定义的两个子程序，并对求解结果进行对比分析。

（3）绘制以上两个方法所求得的方程解的数据分布图。

题目分析：初看题目，MY GOD ！这辈子没见过这么复杂奇葩的方程组。

而脑袋里那些应付考试还可以的少的可怜的Fortran 基础知识，更是直接缴械投降，无地自容了。

不愧是大学，不愧是让无数土木人工科男竞折腰的Fortran 课程设计。

好吧，虽然极度怀疑自己的智商，但不战而屈己之兵又不是我土木人的性格。

打开电脑，摩拳擦掌，Fortran ，老子来了！我们的目的地是，线性方程组的解，通往目的地的道路有好几条，Gauss 大道，矩阵求逆之道，还有两条小路：三角分解与追赶之径。

目前的情况是，小道路黑路不熟，大道倒是有星点的光。

果断走大道嘛。

开发思想嘛，对于解这样一个复杂的线性方程组，聪明的人类是不会傻乎乎自己去做的。

于是，我们把这个繁琐的工作交给任劳任怨的计算机吧，省下我们大好年华去干更多有意义的事情。

数值计算方法FORTRAN程序—函数插值和方程组求解

Difference between the Original Vector and Matrix*P:
3.552713678800501E-015 7.105427357601002E-015 0.000000000000000E+000
1.421085471520200E-014 1.065814103640150E-014 0.000000000000000E+000
FLAG TO CONTROL WHETHER TO PRINT DIFF-CHART:
1
输出数据文件：output2.dat
Nn(x)= 38.500000000000000
Forward difference diagram:
1.000000000
16.00000000 15.00000000
81.00000000 65.00000000 50.00000000
二.
编制解Ax=b的通用子程序。
1)列主元消去法；
流程图（图4、5）
程序（程序四、五）
计算实例
输入数据文件input.dat：
COEFFICIENTS MATRIX IS:
0. 0. 37. 4. 0.
46. 400. 0. 110. 0. 38.
0. 0. 200. 0. 40. 55.
B)程序（程序四、六）
C)计算实例
输入数据文件input.dat：（与上面列主元消去法时相同）
输出数据文件output.dat：
Answer to the simutaneous linear equations is:
-4.312691791318356E-001 -6.207237146208434E-002 1.150000000000000E-001

fortran数值计算基础

数值计算基础目录实验一直接法解线性方程组的 (2)实验二插值方法 (11)实验三数值积分 (5)实验四常微分方程的数值解 (7)实验五迭代法解线性方程组与非线性方程 (9)实验一直接法解线性方程组一、实验目的掌握全选主元消去法与高斯-塞德尔法解线性方程组。

二、实验内容分别写出Guass 列选主元消去法与追赶法的算法，编写程序上机调试出结果，要求所编程序适用于任何一解线性方程组问题，即能解决这一类问题，而不是某一个问题。

实验中以下列数据验证程序的正确性。

1、用Guass 列选主元消去法求解方程组⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--5.58.37.33.47.11.85.16.93.51.53.25.2321x x x2、用追赶法求解方程组⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡-=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡-----000010210000210000210000210000254321x x x x x 三、实验仪器设备与材料主流微型计算机四、实验原理1、Guass 列选主元消去法对于AX =B1）、消元过程：将（A|B ）进行变换为)~|~(B A ，其中A ~是上三角矩阵。

即：⎪⎪⎪⎪⎪⎭⎫⎝⎛→⎪⎪⎪⎪⎪⎭⎫ ⎝⎛n nn nnn nnn n nn b a b a b a a b a a a b a a a b a a a0010122111221222221111211 k 从1到n-1a 、列选主元选取第k 列中绝对值最大元素ik ni k a ≤≤max 作为主元。

b 、换行ik ij kj b b n k j a a ⇔+=⇔,,1,c 、归一化kkk k kj kk kj b a b n k j a a a ⇒+=⇒/,,1,/d 、消元nk i b b a b n k j n k i a a a a i k ik i ij kj ik ij ,,1,,,1;,,1, +=⇒-+=+=⇒-2）、回代过程：由)~|~(B A 解出11,,,x x x n n -。

fortran常见错误及其原因

常见fortran错误1.Incrementally linked image--PC correlation disabled.！编译终止2. forrtl:severe(157):Program Exception - access violationbounds and /warn:argument_checking options set, to see if theproblem is an out-of-bounds memory reference or a argument mismatch that causes data to be treated as anaddress.Other causes of this error include:severe(64):input conversion error, unit 2, file D:\FORTRAN2\testi！文件testi正在读写，直到读写到2时错误。

举例：程序想读写整数，却碰到变量故终止。

4 error LNKZOOI :unresolved external symbol _ SN @ 4 fatal error LNKllZO :1 unresolved externals!出现了未指定的外部函数符号Sn。

这是因为在函数子程序中错把函数名Sn写成了ns。

根据错误信息中的提示，用户在编辑窗口寻找有错位置进行修改。

连接错误往往出现在有函数调用或子程序调用的程序中，常见的错误性质有：未定的函数符号、找不到主程序或子程序、实参与虚参的个数不一致等。

注意：连接错误只给出错误代号和错误性质，不给出具体语句的行号。

5 :: error FOR229O :implicit type for 1 detected between 1 and = C :\ abc．: error FOR33Og :undefined label 10编译系统提示用户：在程序的第5行，变量i未经类型说明；在程序的第H行，标号10未定义。

病态线性方程组的求解

《科学与工程计算》实验报告学号：姓名：一、实验内容：考虑方程组Hx=b 的求解，其中系数矩阵H 为Hilbert 矩阵，nj i j i h h H j i n n j i ,,2,1,,11,)(,, =-+==⨯这是一个著名的病态问题。

通过首先给定解（例如取为各个分量均为1）再计算出右端b 的办法给出确定的问题。

实验要求：（1）选择问题的维数为6，分别用Jacobi 迭代法、GS 迭代法和SOR 迭代法求解方程组，其各自的结果如何？将计算结果与问题的解比较，结论如何？（2）逐步增大问题的维数，仍然用上述的方法来解它们，计算的结果如何？计算的结果说明了什么？（3）讨论病态问题求解的算法。

二、程序设计的基本思想、原理和算法描述：1、算法 Jacobi 迭代法若A 为稀疏矩阵，只需遍历非零元素GS 迭代法若A 为稀疏矩阵，只需遍历非零元素每步迭代计算量相当于一次矩阵与向量的乘法；不需要保留上一步的迭代解，与Jacobi 迭代法计算量一样。

SOR 迭代法（稠密矩阵）2、函数组成double max(double array[100]) 求数组中的最大值函数3、输入/输出设计对于方程组Hx=b 的求解，系数矩阵H 为Hilbert 矩阵，矩阵中的数由下列函数生成。

nj i j i h h H j i n n j i ,,2,1,,11,)(,, =-+==⨯X*取一个特解[1,1,1, (1)b 数组由矩阵的每行元素相加得来。

4、符号名说明double c[100] 用来存储第k+1次和第k 次迭代结果的差值的绝对值 double x[100] 第k+1次迭代结果，储存解数组 double x0[100] 初始向量double r 第k+1次和第k 次迭代结果的差值的绝对值的最大值 double sum 矩阵方程变换后右侧值的和 int k 迭代次数double a[100][100] 存储Hilbert 矩阵 double b[100] 存储b 向量三、源程序及注释：Jacobi 迭代法#include<iostream> #include<math.h> #include <iomanip> #include <stdio.h> using namespace std; int n;double max(double array[100])//求最大值函数 {double a=array[0]; int i;for(i=1; i<n; i++) {if(a<array[i]) a=array[i]; return a; } }double c[100]= {0.0};double x[100]= {0.0}; //第k+1次迭代结果，储存解数组 double x0[100]= {0.0}; //初始向量double r,sum=0; int main(){double s=0;int i,k,j;//k为迭代次数double a[100][100];double b[100]= {0.0};cout<<"请输入维数："<<endl;cin>>n;cout<<"输出a数组："<<endl;double max(double array[100]);for(i=0; i<n; i++){for(j=0; j<n; j++){a[i][j]=1.0/(i+j+1);printf("%10.6lf ",a[i][j]);//矩阵中的数精确到六位}cout<<endl;}cout<<"输出b数组："<<endl;for(i=0; i<n; i++){for(j=0; j<n; j++){b[i]+=a[i][j];//矩阵每行的和}cout<<b[i]<<" ";}cout<<endl;cout<<"输入精度："<<endl;cin>>s;for(k=1;; k++){for(i=0; i<n; i++){for(j=0; j<n; j++){sum=a[i][j]*x0[j]+sum;}x[i]=x0[i]+((b[i]-sum)/a[i][i]);//雅克比迭代方法计算方式c[i]=fabs(x[i]-x0[i]);//求差值的绝对值x0[i]=x[i];sum=0;r=max(c);if(r<s)//输出迭代次数{for(i=0; i<n; i++)cout<<"x"<<i<<" = "<<x[i]<<setprecision(10)<<endl;cout<<"迭代次数："<<k<<endl;return 0;}}return 0;}GS迭代法#include<iostream>#include<math.h>#include <iomanip>#include <stdio.h>using namespace std;int n;double max(double array[100])//求最大值函数{double a=array[0];int i;for(i=1;i<n;i++){if(a<array[i])a=array[i];}return a;}main() {double s=0;double max(double array[100]);double c[100]={0.0};double x[100]={0.0};//第k+1次迭代结果，储存解数组double x0[100]={0.0};//初始向量int i,k,j;double r,sum=0;double a[100][100];double b[100]= {0.0};cout<<"请输入维数："<<endl;cout<<"输出a数组："<<endl;for(i=0; i<n; i++){for(j=0; j<n; j++){a[i][j]=1.0/(i+j+1);printf("%10.6lf ",a[i][j]);//矩阵中的数精确到六位}cout<<endl;}cout<<"输出b数组："<<endl;for(i=0; i<n; i++){for(j=0; j<n; j++){b[i]+=a[i][j];//矩阵每行的和}cout<<b[i]<<" ";}cout<<endl;cout<<"输入精度："<<endl;cin>>s;for(k=1;;k++){for(i=0;i<n;i++){for(j=0;j<n;j++){sum=a[i][j]*x0[j]+sum;}x[i]=x0[i]+((b[i]-sum)/a[i][i]);//gs迭代方法计算方式c[i]=fabs(x[i]-x0[i]);//求差值的绝对值x0[i]=x[i];sum=0;}r=max(c);if(r<s)//输出迭代次数{for(i=0;i<n;i++)cout<<"x"<<i<<" = "<<x[i]<<setprecision(10)<<endl;cout<<"迭代次数："<<k<<endl;return 0;}}}SOR迭代法#include<iostream>#include<math.h>#include <iomanip>#include <stdio.h>using namespace std;int n;double max(double array[100]){double a=array[0];int i;for(i=1;i<n;i++){if(a<array[i])a=array[i];}return a;}int main() {double s=0,w=0;int i,k,j;//k为迭代次数double max(double array[100]);double c[100]= {0.0}; //double x[100]= {0.0}; //第k+1次迭代结果，储存解数组double x0[100]= {0.0}; //初始向量int i,k,j;double r,sum=0;double a[100][100];double b[100]= {0.0};cout<<"请输入维数："<<endl;cin>>n;cout<<"输出a数组："<<endl;for(i=0; i<n; i++){for(j=0; j<n; j++){a[i][j]=1.0/(i+j+1);printf("%10.6lf ",a[i][j]);//cout<<a[i][j]<<" ";}cout<<endl;}cout<<"输出b数组："<<endl;for(i=0; i<n; i++){for(j=0; j<n; j++){b[i]+=a[i][j];}cout<<b[i]<<" ";}cout<<endl;cout<<"输入精度："<<endl;cin>>s;cout<<"输入松弛因子："<<endl;cin>>w;for(k=1;;k++){for(i=0;i<n;i++){for(j=0;j<n;j++){sum=a[i][j]*x0[j]+sum;}x[i]=x0[i]+(w*(b[i]-sum)/a[i][i]);//sor迭代方法的计算公式c[i]=fabs(x[i]-x0[i]);//求差值的绝对值x0[i]=x[i];sum=0;}r=max(c);if(r<s)//输出迭代次数{for(i=0;i<n;i++)cout<<"x"<<i<<" = "<<x[i]<<setprecision(10)<<endl; cout<<"迭代次数："<<k<<endl;return 0;}}}四、运行输出结果：JacobiGSSOR五、结果比较分析：说明：Hx=b，H矩阵可以由函数hi,j=1/(i+j-1)直接由程序生成，为了设定参考解，我们先设x为分量全1的向量，求出b,然后将H和b作为已知量，求x，与设定的参考解对比。

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

PROGRAM TEXT2
PARAMETER(N=3)
REAL,DIMENSION(N,N)::A,B
REAL,DIMENSION(N)::P,Q,E,F
REAL NORM1,NORM2,COND
OPEN(1,FILE='XISHU.TXT')
READ(1,*)((A(I,J),J=1,N),I=1,N)
READ(1,*)(P(I),I=1,N)
READ(1,*)(Q(I),I=1,N)
CLOSE(1)
B=A
CALL triangle(A,P,E,N)
PRINT *,"向量元素b3为0.52时，方程组的解是:"
WRITE(*,100)(E(I),I=1,N)
CALL triangle(A,Q,F,N)
PRINT *,"向量元素b3为0.53时，方程组的解是:"
WRITE(*,100)(F(I),I=1,N)
100 FORMA T(1X,F6.1)
CALL FANSHU(B,N,NORM1)
CALL GAUSSJ(B,N)
CALL FANSHU(B,N,NORM2)
COND=NORM1*NORM2
PRINT *,"矩阵A的条件数为:"
WRITE(*,*)COND
END
SUBROUTINE GAUSSJ(A,N)
INTEGER N,NMAX
REAL A(N,N),B(N)
PARAMETER (NMAX=50)
INTEGER I,ICOL,IROW,J,K,L,L1,INDXC(NMAX),INDXR(NMAX),IPIV(NMAX) REAL BIG,DUM,PIVINV
DO J=1,N
IPIV(J)=0
ENDDO
DO I=1,N
BIG=0
DO J=1,N
IF(IPIV(J)/=1)THEN
DO K=1,N
IF(IPIV(K)==0)THEN
IF(ABS(A(J,K))>=BIG)THEN
BIG=ABS(A(J,K))
IROW=J
ICOL=K
ENDIF
ELSE
IF(IPIV(K)>1)THEN
PAUSE 'SINGULAR MA TRIX IN GAUSSJ'
ENDIF
ENDIF
ENDDO
ENDIF
ENDDO
IPIV(ICOL)=IPIV(ICOL)+1
IF(IROW/=ICOL)THEN
DO L=1,N
DUM=A(IROW,L)
A(IROW,L)=A(ICOL,L)
A(ICOL,L)=DUM
ENDDO
DUM=B(IROW)
B(IROW)=B(ICOL)
B(ICOL)=DUM
ENDIF
INDXR(I)=IROW
INDXC(I)=ICOL
IF(A(ICOL,ICOL)==0.)PAUSE 'SINGULAR MA TRIX IN GAUSSJ' PIVINV=1./A(ICOL,ICOL)
A(ICOL,ICOL)=1
DO L=1,N
A(ICOL,L)=A(ICOL,L)*PIVINV
ENDDO
B(ICOL)=B(ICOL)*PIVINV
DO L1=1,N
IF(L1/=ICOL)THEN
DUM=A(L1,ICOL)
A(L1,ICOL)=0
DO L=1,N
A(L1,L)=A(L1,L)-A(ICOL,L)*DUM
ENDDO
B(L1)=B(L1)-B(ICOL)*DUM
ENDIF
ENDDO
DO L=N,1,-1
IF(INDXR(L)/=INDXC(L))THEN
DO K=1,N
DUM=A(K,INDXR(L))
A(K,INDXR(L))=A(K,INDXC(L))
A(K,INDXC(L))=DUM
ENDDO
ENDIF
ENDDO
ENDDO
END
SUBROUTINE FANSHU(B,N,AM) DIMENSION B(N,N),D(N)
REAL B,D,AM
DO J=1,N
D(J)=0.0
DO I=1,N
D(J)=D(J)+ABS(B(I,J))
ENDDO
ENDDO
AM=0.0
DO J=1,N
IF(AM<D(J)) AM=D(J)
ENDDO
END
SUBROUTINE triangle(A,M,E,N) INTEGER S,P
REAL,DIMENSION(N,N)::A,B,C
REAL,DIMENSION(N)::M,D,E
DO I=1,N
C(1,I)=A(1,I)
B(I+1,1)=A(I+1,1)/C(1,1)
ENDDO
DO I=2,N
DO J=I,N
C(I,J)=A(I,J)
DO S=1,I-1
C(I,J)=C(I,J)-B(I,S)*C(S,J)
ENDDO
ENDDO
J=I
DO P=J,N
B(P+1,J)=A(P+1,J)/C(J,J)
DO S=1,J-1
B(P+1,J)=B(P+1,J)-B(P+1,S)*C(S,J)/C(J,J) ENDDO
ENDDO
ENDDO
DO I=2,N
DO J=1,I-1
C(I,J)=0
ENDDO
ENDDO
DO J=1,N
DO I=1,J
IF(I==J)THEN
B(I,J)=1
ELSE
B(I,J)=0
ENDIF
ENDDO
ENDDO
D(1)=M(1)
DO I=2,N
D(I)=M(I)
DO S=1,I-1
D(I)=D(I)-B(I,S)*D(S)
ENDDO
ENDDO
E(N)=D(N)/C(N,N)
DO I=N-1,1,-1
E(I)=D(I)/C(I,I)
DO S=I+1,N
E(I)=E(I)-C(I,S)*E(S)/C(I,I)
ENDDO
ENDDO
END。