数值线性代数课程设计—超定方程组的求解

第一篇：数值线性代数课程设计—超定方程组的求解

《数值线性代数课程设计》

专业： 信息与计算科学

班级： 13405011 学号： 1340501123 姓名： 实验日期：报告日期：实验地点：邢耀光 数理学院五楼机房

2016.05.09

2015.05.13

超定方程组的求解

邢耀光

（班级：13405011 学号1340501123）

摘要：在实验数据处理和曲线拟合问题中，求解超定方程组非常普遍。比较常用的方法是最小二乘法。形象的说，就是在无法完全满足给定条件的情况下，求一个最接近的解。最小二乘法（又称最小平方法）是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。

关键字：最小二乘问题，残量，超定方程组，正则化方程组，Cholesky分解定理。

正文：

最小二乘法的背景：

最小二乘法（又称最小平方法）是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。最小二乘法还可用于曲线拟合。其他一些优化问题也可通过最小化能量或最大化熵用最小二乘法来表达。最小二乘法经常运用在交通运输学中。

交通发生预测的目的是建立分区产生的交通量与分区土地利用、社会经济特征等变量之间的定量关系，推算规划年各分区所产生的交通量。因为一次出行有两个端点，所以我们要分别分析一个区生成的交通和吸引的交通。

最小二乘问题：

最小二乘问题多产生于数据拟合问题。例如，假定给出m个点t1,...,tm和这m个点上的实验或观测数并假定给出在ti上取值的n个已知函数1(t),...,n(t)。考虑i 的线性组合，f(x;t)x11(t)x22(t)...xnn(t)，（1）

我们希望在t1,...,tm点上f(x;t)能最佳的逼近y1,...,ym这些数据。为此，若定义残量 据y1,...,ymj1

则问题成为：估计参数x1,...,xn，使残量r1,...,rm尽可能地小。（2）式可用矩阵-向量形式表示为

ri(x)yixjj(ti)，i1,...,m，（2）

n r(x)bAx，（3）其中

1(t1)n(t1)y1A, b，(t)(t)y1mnmm

TT)r(x(x,...x，x)(r(x),...,r(x)).1nm1

当mn时，我们可以要求r(x)0，则估计x的问题就可以用第一章中讨论的方法解决。当mn时，一般不可能使所有残量为零，但我们可要求残向量r(x)在某种范数意义下最小。最小二乘问题就是求x使残向量r(x)在2范数意义下最小。

定义1：给定矩阵ARmn及向量bRm，确定xRn，使得

bAx2r(x)2minr(y)2minAyb2.（4）

yRnyRn这就是所谓的最小二乘问题，简称为LS问题，其中的r(x)常常被称为残向量。

在所讨论的最小二乘问题中，若r线性依赖于x，则称其为线性最小二乘问题：若r非线性依赖于x，则称其为非线性最小二乘问题。

最小二乘问题的解x又可称做线性方程组

Axb,ARmn

（5）的最小二乘解，即x在残向量r(x)bAx的2范数最小的意义下满足方程组（5）。当mn时称（5）式为超定方程组。

定理1:（Cholesky分解定理）若ARnn对称正定，则存在一个对角元均为正数的下三角阵LRnn，使得

ALL.（6）（6）式称为Cholesky分解，其中的L称作A的Cholesky因子。

因此，若线性方程组Axb的系数矩阵是对称正定的，则我们自然可按如下的步骤求其解：

（1）计算A的Cholesky分解：ALL ；

（2）求解Lyb得y ；

（3）求解Lxy得x； 简单而实用的方法是直接比较ALL两边的对应元素来计算L。设

TTTTl11ll2122.Llllnnn1n2T比较ALL两边对应的元素，得关系式

aij 首先，由a11l11，得

l11再由ai1l11li1，得

li1ai1l11，i1,...,n.这样便得到了矩阵L的第一列元素。假定已经算出L的前k1列元素，由

akk得 2lp1jipjpl，1jin（7）

a11.lp1k2kp，122 lkkakklkp.（8）

p1k1再由

得

aikliplkpliklkk，ik1,...,n，p1k1k1 likaikliplkplkk，ik1,...,n.（9）

p1这样便求出了L的第k列元素。这种方法称为平方根法。

记最小二乘解的解集为LS，即

定理 ATAxATb.（10）

方程组（10）常常被称为最小二乘问题的正则化方程组或法方程组，它是一个含有n个变量和n个方程的线性方程组。在A的列向量线性无关的条件下，AA对称正定，故可用平方根法求解方程组（6），这样，我们就得到了求解最小二乘问题最古老的算法———正则化方法，其基本步骤如下：

（1）计算CAA, dAb;

（2）用平方根法计算C的Cholesky分解：CLL;（3）求解三角方程组Lyd和Lxy.TTTTT

2：xLS 当且仅当

LSxRn:x是LS问题(3)的解，实验 ：

一：超定方程组的求解

原理：设A是mn阶矩阵mn，则线性方程组Axb为超定方程组，这里xR,bR。如

mm果A的秩为n，则称A为列满秩矩阵。超定方程组的解满足法方程AAxAb，该解使得

TTbAx 22min，称之为最小二乘解。

11 题目： 111TT1.11.12121.21.2221.31.3x3

21.41.441.51.525

用正则化方法求解，要求:（1）BLL 不得使用MathCAD指令Cholesky；（2）BLL使用MathCAD指令Cholesky。

11解：（1）A111 1.11.126.58.5551.21.22T8.5511.375 1.31.32 , 则 BAA6.51.41.428.5511.37515.2981.51.5215B21,20.5L2.907 , gATbLB2.236, 211111L1128.25

L31 B3123.824 , L22B22L210.316 , L11 4 LB32L31L2132L0.822 , LL2233B3331L320.037 , 22

2.236002.2362.9073.824 即 L2.9070.31600.037 , LT00.3160.822, 3.8240.822000.037

6.70810 则yL1g3.162xLT1y10 ,  , 9.27310132.4781011

x即为所求的最小二乘解。

11.11.1211.21.22（2）A2.2360011.31.32cholesky(B)2.9070.316011.41.423.8240.8220.03711.51.52，2.236002.2362.9073.824 则 L2.9070.3160，LT00.3160.8223.8240.8220.037000.037，6.70810 则yL1g3.162xLT1y10 , 

9.27310132.4781011,x即为所求的最小二乘解。

二：已知如下数据： xi0.00.20.40.60.81.01.2yi0.91.92.83.34.05.76.5 利用最小二乘法拟合曲线 ya1xa2.0.00.90.21.9解：令B0.00.20.40.60.81.01.20.42.80.91.92.83.34.05.76.5 ,x0.6，y3.3 0.81.04.05.71.26.510.010.210.41 则A10.6XATAATy0.8434.571，即p(x)0.8434.571x， , 10.811.011.2故最小二乘法拟合曲线为y4.571x0.843.程序附录： 一；

11.11.12111.21.22256.58.55A11.31.32b3AxbBATA gATb, B6.58.5511.37511.41.428.5511.37515.298,45, ，11.51.52 f(B)nrows(B)

Lidentityn()fork1nLkkBkkifk1k1LkkBkkLkp2otherwisep1forik1nBLkikiLifk1kk(break)ifknk1BikLipLkpLikp1LotherwisekkL

15g20.5, , 28.25，002.23602.23602.2362.9073.824Tf(B)2.9070.3160L2.9070.3160L00.3160.82200.0373.8240.8220.037,Lf(B),3.8240.8220.037,0，6.708103.16210yx119.27310132.4781011TyLg , xLy ,, ，二；

120.00.20.40.60.81.01.2TT B , xB , yB, 0.91.92.83.34.05.76.5

00.90.21.90.42.8x0.6 ,y3.3 , x0, x0.2,nrows(x),n120.8415.71.26.507,i1n , Ai11，.111Ax,AXy,A1i2i111p(s)0.8434.571x.0.20.410.843,p(s)XT1

TTAy ,X0.6,XAA4.571s0.811.2

心得体会：

通过本次的课程设计，让我学会了很多，学会了简单的MathCAD 软件的用法。让我更加深刻了解最小二乘问题，和以往对知识的疏忽得以补充。不仅掌握了学习的知识，而其还可以培养和熟练使用资料，运用工具书的能力，把我们所学的课本知识与实践结合起来，起到温故而知新的作用。

参考文献：

1，《数值线性代数》（第二版）北京大学出版社 徐树方，高立，张平方，编著。2013.01

email：974671870qq.com

1340501123：邢耀光

2016.05.13

第二篇：数据结构课程设计-迷宫求解范文

数据结构课程设计

迷宫求解

学院：湖北工业大学计算机学院

教师：沈华老师

班级：12网络工程1班

学号：1210322118

姓名：饶进阳

时间：2013年12月22日

目

录

问题描述

.........................................................思路解析

.........................................................程序流程

.........................................................核心代码

.........................................................源程序代码........................................................程序运行实例.....................................................课设总结

.........................................................3 4 5 6 12 14 迷宫求解

问题描述：

可以输入一个任意大小的迷宫数据，用非递归的方法求出一条走出迷宫的路径，并将路径输出；

要求：

在上交资料中请写明：存储结构、基本算法（可以使用程序流程图）、源程序、测试数据和结果、算法的时间复杂度、另外可以提出算法的改进方法；

比如这是一个迷宫

电脑找出的出路

思

路

设定当前位置的初值为入口位置： do｛

若当前位置可通，则｛

将当前位置插入栈顶；} 若该位置是出口位置，则结束；

否则切换当前位置的东邻块为新的当前位置； 否则｛

若栈不空且栈顶位置还有其他方向未被探索，则设定新的当前位置为沿顺时针方向旋转找到的栈顶位置的下一相邻块；} 若{ 栈不空但栈顶位置的四周均不可通，则删去栈顶位置；} 若{ 栈不空，则重新测试新的栈顶位置，直至找到一个可通的相邻块或出栈至栈空；｝ ｝while(栈不空)｛栈空说明没有路径存在｝

PS：类似动态求解方法

程序流程

第一次用visio做程序框图，所以只把程序的大概流程画出来了。

核心代码

//栈相关操作 int initlStack(Stack *S)int pushStack(Stack S, coordinate e)int popStack(Stack S, coordinate *e)int getTop(Stack S, coordinate *e)void show(Stack S)

//创建一个迷宫 int creatMaze(Maze *m)

//打印出一个迷宫 void showMaze(Maze *m)

//求当前结点的下一个通路

coordinate passNext(Maze *m, int i, int j)//求迷宫路径

int solveMaze(Maze *m, Stack S)//模拟出路径

void showRoad(Maze m, Stack S)

程序源码：

#include #include

#define MAX 100

#define SIZE sizeof(Node)//**************************** //栈

typedef struct { int x;//行

int y;//列 }coordinate;//坐标

typedef struct node{ coordinate data;struct node *next;}Node, *Stack;

//栈的相关操作

//栈初始化

//带头结点的链栈

int initlStack(Stack *S){(*S)=(Node *)malloc(SIZE);if((*S)== NULL){

return 1;}(*S)->next = NULL;return 0;}

//入栈

int pushStack(Stack S, coordinate e){ Node *p;

p =(Node *)malloc(SIZE);p->data.x = e.x;p->data.y = e.y;p->next = S->next;S->next = p;return 0;}

//出栈

int popStack(Stack S, coordinate *e){ Node *p;if(S->next == NULL){

printf(“nStack is empty!n”);

return 2;} e->x = S->next->data.x;e->y = S->next->data.y;p = S->next;S->next = p->next;free(p);return 0;}

//读取栈顶

int getTop(Stack S, coordinate *e){ e->x = S->next->data.x;e->y = S->next->data.y;return 0;}

//显示

void show(Stack S){ Node *p;p = S->next;while(p!= NULL){

printf(“%d %d <-”, p->data.x, p->data.y);

p = p->next;} printf(“startn”);} //***************************

//*************************** //迷宫

typedef struct { int row;int col;int arr[MAX][MAX];}Maze;

//创建一个迷宫

int creatMaze(Maze *m){ int i, j;printf(“nplease input Maze's row and col:n”);scanf(“%d%d”, &(m->row), &(m->col));printf(“nplease input the Maze's message:(0 is ok),(1 is no)n”);for(i = 0;i <= MAX1;j++){

m->arr[i][j] = 1;

} } for(i = 1;i <= m->row;i++){

for(j = 1;j <= m->col;j++){

scanf(“%d”, &(m->arr[i][j]));

} } return 0;}

//显示迷宫信息

void showMaze(Maze *m){ int i, j;printf(“nthe Maze you hava input:n”);for(i = 1;i <= m->row;i++){

for(j = 1;j <= m->col;j++){

printf(“%d ”, m->arr[i][j]);

}

printf(“n”);}

printf(“nlike this:n”);for(i = 1;i <= m->row;i++){

for(j = 1;j <= m->col;j++){

if(m->arr[i][j]!= 0){

printf(“■ ”);

}

else {

printf(“□ ”);

}

}

printf(“n”);} } //求当前结点的下个通路 //顺时针找

coordinate passNext(Maze *m, int i, int j){ coordinate k;if(m->arr[i][j+1] == 0){//东

k.x = i;

k.y = j+1;

return k;}

if(m->arr[i+1][j] == 0){//南

k.x = i+1;

k.y = j;

return k;} if(m->arr[i][j-1] == 0){//西

k.x = i;

k.y = j-1;

return k;} if(m->arr[i-1][j] == 0){//北

k.x = i-1;

k.y = j;

return k;} else {//不存在 k.x =-1;

k.y =-1;

return k;} }

//求解迷宫

int solveMaze(Maze *m, Stack S){ int i, j;int boolean;coordinate a;i = 1;j = 1;do {

if(m->arr[i][j] == 0){

a.x = i;

a.y = j;

m->arr[i][j] = 2;

pushStack(S, a);

if(i == m->row && j == m->col){

return 0;

}

}

a = passNext(m, i, j);

if(a.x!=-1){

i = a.x;

j = a.y;

continue;

}

else if(a.x ==-1){

boolean = popStack(S, &a);

if(2 == boolean){

return 0;

}

getTop(S, &a);

i = a.x;

j = a.y;

} }while(S->next!= NULL);return 0;}

//模拟出路径

void showRoad(Maze m, Stack S){ coordinate e;int i, j;if(S->next == NULL){

printf(“nSorry!you can't go out the Maze!n”);} while(S->next!= NULL){

popStack(S, &e);

m.arr[e.x][e.y] =-1;} for(i = 1;i <= m.row;i++){

for(j = 1;j <= m.col;j++){

if(m.arr[i][j] >= 0){

printf(“■ ”);

}

else {

printf(“□ ”);

}

}

printf(“n”);} } //*************************** //主函数 int main(){ Maze m;Stack S;printf(“※ 欢迎玩耍迷宫求解游戏 ※n”);initlStack(&S);creatMaze(&m);showMaze(&m);solveMaze(&m, S);printf(“nthe root is here:n”);show(S);showRoad(m, S);return 0;}

程序运行结果：

第一步：输入迷宫行列大小

第二步：输入迷宫信息

第四步：程序会自动求出路径

课 设 总 结

敬爱的沈华老师您好，感谢您这半年对我们的教诲，说实话，很喜欢上您的课，您上课的风采深深的吸引了我，每次上您的课，我都听得特别认真。好吧，言归正传，谈谈这次课程设计的感想。

首先看到这道题的时候，真心没多少思路，我抱着试一试的态度，去网上搜索相关的算法。网上资源蛮多的，刚开始就搜索到了，严蔚敏老师的相关算法讲解视频，虽然她讲的是思想（就是他并没有源程序的实现），但是足够让我理解了基本算法流程。

其实先开始都不知道怎么下笔，但是我还是想着从基本栈开始写吧，慢慢写着，程序大概思想框架就成型了，最后在实现求解迷宫路径的时候，着手用笔在草稿纸上面模拟，然后才敲到编译器中，最好几经调试，终于得到了想要的结果。

做完这次课设后，深刻的体会到一些道理，在信息时代，解决问题的最好方法是运用现代的信息资源。并且在实际中，开始没思路的事，慢慢完成小部分，那么你的总体思维会渐渐地成型。

做事做人，平心气和的干，总会得到成功！

第三篇：数值线性代数课设课件资料

数值线性代数课程设计报告

姓名:陶英 学号:081410124

任课教师:杨熙

南京航空航天大学

2016 年 6 月 22日

求解线性方程组的三种迭代法及其结果比较

摘要

当今的环境下，数值计算越来越依赖于计算机。大规模科学计算和工程技术中许多问题的解决，最终归结为大型稀疏线性方程组的求解，其求解时间在整个问题求解时间中占有很大的比重，有的甚至达到80%。由于现今科学研究和大型项目中各种复杂的可以对计算精度和计算速度的要求越来越高。因此，作为大规模科学计算基础的线性代数方程组的高效数值求解引起了人们的普遍关注。这种方程组的求解一般采用迭代法。

关于迭代法，是有很多种解决公式的：Jacobi，G-S和超松弛迭代法。这三种方法的原理大致相同，Jacobi需要给定初向量，G-S则需要给定初值，超松弛法是对Guass-Seidel迭代法的加权平均改造。而本文则是对大型稀疏线性方程组迭代求解与三种迭代法（Jacobi，Gauss-Seidel和超松弛迭代法）的收敛速度与精确解的误差比较做出研究。

关键词：Jacobi迭代法；Gauss-Seidel迭代法；SOR迭代法；线性方程组 方法与理论的叙述

1.1迭代法简介

1.Jacobi迭代法：

对于非奇异线性方程组Ax=b，令A=D-L-U，其中

则原方程组可改写为：

其中

给定初始向量：

由（2.2）可以构造迭代公式：

其分量形式为：

2.2）（2.Guass-Seidel迭代法： 类似于Jacobi迭代法，给定初值：

令

则得到Guass-Seidel公式：

其分量形式为：

3.超松弛迭代法（SOR 迭代法）：

SOR迭代法是对Guass-Seidel迭代法的加权平均改造，即

为Guass-Seidel迭代解，即

它的分量形式为：

其中ω称为松弛因子，当ω>1时称为超松弛；当ω<1时叫低松弛；ω=1时就是

Guass-Seidel迭代。

上述三种经典迭代法收敛的充分必要条件是迭代矩阵谱半径小于1。

谱半径不易求解，而在一定条件下，通过系数矩阵A的性质可判断迭代法的收敛性。定理1：

若系数矩阵A是严格对角占优或不可约对角占优，则Jacobi迭代法和Gauss-Seidel迭代法均收敛。定理2：

（1）SOR迭代法收敛的必要条件是0

（2）若系数矩阵A严格对角占优或不可约对角占优且0

2.1问题

考虑两点边值问题：

d2ydya,0a12 dxdxy(0)0,y(1)1.1a(1e)ax 容易知道它的精确解为：y1/1ex为了将微分方程离散，把[0,1]区间n等分，令h=1/n，xiih,i1,2,...n1，得到差分方程

(h)yi1(2h)yiyi1ah2，从而得到迭代方程组的系数矩阵A。

对=1，a=1/2,n=100,分别用jacobi，G-S，超松弛迭代法分别求线性方程组的解，要求4位有效数字，然后比较与精确解的误差。

对=0.1，=0.01，=0.001，考虑同样问题。

1.方程的表示及存储

由于本题中线性方程组的系数矩阵为三对角矩阵，所以可以采用紧缩方法存储，即

然后在矩阵乘法时对下标处理一下即可。但是考虑到三种迭代方法的一般性，且本题中n=200并不是很大，所以实验中并没有采用紧缩存储，而是采用了直接存储。2.边值条件的处理

由于差分得到的方程组的第一行和最后一行中分别出现了边值y(0)与y(1)作为常数项，因此要在常向量的第一项和最后一项作一些修改：

3.迭代终止条件

首先确定要求的精度tol，我们希望当

则停止迭代。对于迭代格式，若

且，则迭代序列的

第k 次近似解和精确解之间有估计式由题目要求知我们需要有

。，而由上面的迭代估计，只要，即取为，因此最后令迭代终止条件为

即可。而本题中q可近似

4.SOR 迭代中最佳松弛因子的选取

由于SOR 迭代法的效果和其松弛因子w的选取有关，所以有必要选取合适的松弛因子。当选择最佳松弛因子

时，SOR 方法的迭代速度最快。

Matlab实现：

迭代矩阵是n-1阶的，不是n阶；

等号右端向量b的最后一项，不是ah^2，而是ah^2-eps-h

2.2精确解

1ay(1e)ax 1/1ex带入a=1/2,=1 代码： >> clear >> x=linspace(0,1);truy=(1-0.5)/(1-exp(-1/1))*(1-exp(-x./1))+x.*0.5;figure;plot(x,truy,'g','LineWidth',1.5);hold on;Grid

图：

2.3三种迭代法

Jacobi法：代码见附录 Eps=1 结果：

迭代次数k：22273 结果与精确解的比较图（绿色粗线是精确解，黑色细线是迭代结果）

Eps=0.1 结果：

迭代次数k：8753 结果与精确解的比较图（绿色粗线是精确解，黑色细线是迭代结果）

Eps=0.01 结果：

迭代次数k：661 结果与精确解的比较图（绿色粗线是精确解，黑色细线是迭代结果）

G-S迭代法：代码见附录 Eps=1 结果：

迭代次数k：11125 结果与精确解的比较图（绿色粗线是精确解，黑色细线是迭代结果）

Eps=0.1 结果：

迭代次数k：4394 结果与精确解的比较图（绿色粗线是精确解，黑色细线是迭代结果）

Eps=0.01 结果：

迭代次数k：379 结果与精确解的比较图（绿色粗线是精确解，黑色细线是迭代结果）

超松弛法：代码见附录 Eps=1 w=1.56 结果：

迭代次数k：3503 结果与精确解的比较图（绿色粗线是精确解，黑色细线是迭代结果）

Eps=0.1 w=1.56 结果：

迭代次数k：1369 结果与精确解的比较图（绿色粗线是精确解，黑色细线是迭代结果）

Eps=0.01 w=1.56 结果：

迭代次数k：131 结果与精确解的比较图（绿色粗线是精确解，黑色细线是迭代结果）分析讨论及心得体会

3.1三种方法的比较

Jacobi、G-S、超松弛法，三者都能够取得对精确解的良好逼近，但是，在相同的精度条件下，三者的收敛速度是不一样的，jacobi

3.2心得体会

这次课程设计，平时感觉挺简单的那些枯燥单调的代码和数学公式，真正到了自己运用的时候却无从下手，但是，解决问题的过程恰是不断学习的过程：数学算法转换为代码的过程要对题目有深入的了解，然后对程序函数定义还要有一定的掌握能力，通过这个的过程让我巩固了自己的数学知识，对数学专业知识和MATLAB的操作有了更深的体会。

课程设计中遇到的问题只凭自己苦思冥想是不能全部解决的，这是同学老师的建议和网络给了我很大的帮助。遇到自己解决不了的问题时，多多向老师同学请教，或许问题就能迎刃而解。

4参考文献

[1]徐树方.数值线性代数.北京:北京大学出版社,1995.[2]马昌凤.现代数值分析.北京:国防工业出版社.2013.[3]刘春凤,米翠兰.实用数值分析教程.北京冶金工业出版社.2006

5附录

源代码

1.Jacobi：

function [y,k]=jacobi2(a,eps,h,delta)n=1.0/h;A=ones(n-1);y=zeros(n-1,1);z=zeros(n-1,1);k=0;for i=1:n-1 for j=1:n-1 A(i,j)=0;end end

for i=1:n-1 A(i,i)=-(2*eps+h);end

for i=1:n-1 for j=1:n-1 if i==j+1 A(i,j)=eps;end

if i==j-1 A(i,j)=eps+h;end

end end

b=zeros(n-1,1);for i=1:n-2 b(i,1)=a*h^2;end

b(n-1,1)=a*h^2-eps-h;D=zeros(n-1);for i=1:n-1 D(i,i)=A(i,i);end

L=zeros(n-1);for i=1:n-1 for j=1:n-1 if i>j L(i,j)=-A(i,j);end

end end

U=zeros(n-1);for i=1:n-1 for j=1:n-1 if i

end end

B=D(L+U);g=Db;while 1 z=B*y+g;if norm(z-y,inf)

y=z;k=k+1;end

x=linspace(0,1);truy=(1-a)/(1-exp(-1/eps))*(1-exp(-x./eps))+x.*a;figure;plot(100*x,truy,'g','LineWidth',5);hold on;grid hold on;plot(y,'b')

2.G-S: function [y,k]=gs2(a,eps,h,delta)n=1.0/h;A=ones(n-1);y=zeros(n-1,1);z=zeros(n-1,1);k=0;for i=1:n-1 for j=1:n-1 A(i,j)=0;end end

for i=1:n-1 A(i,i)=-(2*eps+h);end

for i=1:n-1 for j=1:n-1 if i==j+1 A(i,j)=eps;end

if i==j-1 A(i,j)=eps+h;end

end end

b=zeros(n-1,1);for i=1:n-2 b(i,1)=a*h^2;end

b(n-1,1)=a*h^2-eps-h;D=zeros(n-1);for i=1:n-1 D(i,i)=A(i,i);end

L=zeros(n-1);for i=1:n-1 for j=1:n-1 if i>j L(i,j)=-A(i,j);end

end end

U=zeros(n-1);for i=1:n-1 for j=1:n-1 if i

end end

B=D(L+U);g=Db;while 1 z=(D-L)U*y+(D-L)b;if norm(z-y,inf)

y=z;k=k+1;end

x=linspace(0,1);truy=(1-a)/(1-exp(-1/eps))*(1-exp(-x./eps))+x.*a;figure;plot(100*x,truy,'g','LineWidth',5);hold on;grid hold on;plot(y,'b')

3.SOR：

function [y,k]=sor(a,eps,h,delta,w)n=1.0/h;A=ones(n-1);y=zeros(n-1,1);z=zeros(n-1,1);k=0;for i=1:n-1 for j=1:n-1 A(i,j)=0;end end

for i=1:n-1 A(i,i)=-(2*eps+h);end

for i=1:n-1 for j=1:n-1 if i==j+1 A(i,j)=eps;end if i==j-1 A(i,j)=eps+h;end

end end

b=zeros(n-1,1);for i=1:n-2 b(i,1)=a*h^2;end

b(n-1,1)=a*h^2-eps-h;D=zeros(n-1);for i=1:n-1 D(i,i)=A(i,i);end

L=zeros(n-1);for i=1:n-1 for j=1:n-1 if i>j L(i,j)=-A(i,j);end

end end

U=zeros(n-1);for i=1:n-1 for j=1:n-1 if i

end end

B=D(L+U);g=Db;Lw=((D-w*L)^-1)*((1-w)*D+w*U);while 1 z=Lw*y+w*(D-w*L)^-1*b;if norm(z-y,inf)

y=z;k=k+1;end

x=linspace(0,1);truy=(1-a)/(1-exp(-1/eps))*(1-exp(-x./eps))+x.*a;figure;plot(100*x,truy,'g','LineWidth',5);hold on;grid hold on;plot(y,'b')

第四篇：无网格数值求解方法学习小结

无网格数值求解方法

——学习小结

一、无网格法的介绍

有限元法存在的那些问题都来源于网格，在用有限元方法处理诸如金属冲压成型、高速冲击、动态裂纹扩展、流固耦合等涉及大变形和移动边界的问题时，由于网格可能发生严重扭曲，往往需要网格重构，不但精度受到了严重影响，计算也大幅度提高，因此有限元方法在这些领域的应用遇到了困难。

直接在有限元基础上对其进行改进，效果自然不会达到最好，于是研究者把革命的对象锁定在了网格上。几经尝试以后，一种基于点集的插值方法被研究者广泛采用，现今的无网格方法，一般就指的是这一类基于点集的数值方法。

无网格方法的位移函数是在点的领域内构造的，并且这些区域是可以重叠的，因此在处理大变形和移动边界等问题时，没有网格的初始划分和重构问题，这不仅有利于这类问题计算精度的提高，还可以减少数值计算难度。

目前已存在十余种无网格方法，它们之间的区别主要在于试函数的选择和微分方程的等效形式。虽然无网格方法对于大变形和移动边界问题具有优势，但其存在收敛性、数值稳定性和效率等问题，因此无网格方法还只能作为有限元方法的补充。

无网格方法基本思想是将有限元法中的网格结构去除，完全代之以一系列的结点排列。

二、求解方法方法

基于位移最小二乘(MLS)近似方法—EFG(Element-free Galerkin Method, Belytschko, 1994)。EFG方法计算稳定，精度较高，是无网格方法中较为成熟的一种 方法。

无网格法就目前来说，仍没有有限元法发展得那么快。而且，大规模地使用无网格法将大大增加计算时间。因此通常只需要在那些不连续、大变形或应力集中区域使用无网格法进行离散，如冲击区域、裂纹扩展区域、大变形区域等，其余区域仍然可采用其他数值方法。

微分方程组

A

(u)=A1(u)A2(u)0 在 内  ...

边界条件

B1(u)B

(u) B(u)20 在上...

等效积分形式 U TAudVTB ud0(*)等效积分弱形式

CTUDudETVFud0

2.1加权余量法

求解域Ω中，若场函数是精确解，则在域Ω中任一点都满足微分方程，同时在边界上任一点都满足边界条件式，此时等效积分形式或等效积分弱形式必

(**)

然严格地得到满足。但是对于复杂的实际问题，这样的精确解往往是很难找到的，因此, 人们需要设法找到具有一定精度的近似解。设u是一个近似解，即为试函数，它可以表示成为一组已知函数或Ritz基函数i的线性组合，即

uxiTaiTa

i1n式中ai为待定系数或Ritz基坐标。

将权函数代入加权余量积分式，由于系数bj的任意性，有

TTTTRadRjAjBad0,j1,2,,m

上式给出了m个方程。用于求解n个待定系数ai。如果mn，则上式是超定的，需要借助于最小二乘法解。对上式进行分部积分得到等效积分弱形式的近似形式

CDadEFad0

TjTTjT2.2伽辽金法

按照对权函数的不同选择就得到不同的加权余量的计算方法并赋以不同的名称。如果取权函数与试函数相同，则称为Galerkin方法。

nTnTTiaidjRBiaid0 RAi1i1Tj我们将会看到，在很多情况下，采用伽辽金法得到的求解方程的系数矩阵是对称的，这是在用加权余量法建立有限元格式时几乎毫不例外地采用伽辽金法的主要原因，而且当存在相应的泛函数时，伽辽金法与变分法往往导致同样的结果。

2.3移动最小二乘近似

构造方法：考虑求解域，其中共有N个结点xi(i1,2,,N)，在各个结

点处有u0(xi)ui，但u(xi)ui。考虑计算点x（对于无网格配点法为结点；对于伽辽金无网格方法为高斯积分点），其邻域x内的近似函数可以写为

)pi(x)ai(x)pT(x)a(x)u(x,xi1m

[x)为Rits基函数，ai(x)为Rits基坐标或待求系数，x式中：pi(x计算点x邻域x内任意点的坐标，它包括x，m是基函数的个数。而

yz]是)[p1(x)pT(x)pm(x)]p2(x，aT(x)[a1(x)a2(x)am(x)]

值得注意的是，在经典Ritz方法中，Ritz基坐标是常数，并且基函数要满足位移边界条件。在式(1)中，基函数要满足如下条件：

)1p1(x )Cn()pi(x

式中：i1,2,,m，Cn()表示在域内具有直到n阶连续导数的函数空间。2.4边界条件

无网格方法的结点形函数多数都不满足关系Njxiij，因此位移边界条件的处理是比较困难的。若采用紧支径向基函数来构造形函数，则可以像一般有限元方法那样来处理位移边界条件。在MLS近似中，若选奇异函数为权函数，则近似函数具有插值特性即Njxiij，因此可以直接施加本质边界条件。对与其他情况，可以借助拉格朗日乘子方法来处理边界条件。

拉格朗日乘子法包括两种，一种是利用边界积分中直接引入边界条件，即

εTσuTfduTpduuTλ+λTu-ud0

三、具体算例

左端固定的悬臂梁，右端面受抛物线剪切载荷作用

主程序：

tic clear;Lx = 20;Ly = 10;young = 210;nu=0.3;q =-1;a = 0;nx = 30;ny = 20;ndivl=10;ndivw=6;dmax=2.89;Dmat =(young/(1-nu^2))*[1 nu 0;nu 1 0;0 0(1-nu)/2];[x,numnod,dm] = mesh1(Lx,Ly,nx,ny,dmax);figure hold on plot(x(1,1:(ny+1)),x(2,1:(ny+1)),'k-','linewidth',3);axis equal;plot(x(1,(ny+1):(ny+1):numnod),x(2,(ny+1):(ny+1):numnod),'k-','linewidth',2);plot(x(1,numnod:-1:(numnod-ny)),x(2,numnod:-1:(numnod-ny)),'k-','linewidth',2);plot(x(1,1:(ny+1):(numnod-ny)),x(2,1:(ny+1):(numnod-ny)),'k-','linewidth',2);%plot(x(1,:),x(2,:),'k.');___axis off;plot(x(1,:),x(2,:),'k.');axis equal;axis off;hold off [xc,conn,numcell,numq] = mesh2(Lx,Ly,ndivl,ndivw);

[nnu,nnt,numT1,numT2] = mesh3(numq,xc,Lx,Ly,a);% nnu---% nnt---% numT1--% numT2--% numq-----quado = 4;[gauss] = gauss2(quado);numq2 = numcell*quado^2;gs = zeros(4,numq2);[gs] = egauss(xc,conn,gauss,numcell);[k]=kjuzhen(numnod,gs,x,dm,dmax,Dmat);rfa=400e12;[ka]=kajuzhen(numnod,nnu,numT1,xc,gauss,x,dm,dmax,rfa);K=k+ka;[f] = fjuzhen(numnod,nnt,numT2,xc,gauss,x,dm,dmax,q,Ly);%fa = zeros(2*numnod,1);%[fa] fajuzhen(nu,young,q,numnod,nnu,numT1,xc,gauss,x,dm,dmax,rfa,Ly);[fa] fajuzhen(nu,young,q,numnod,nnu,numT1,xc,gauss,x,dm,dmax,rfa,Lx,Ly)F=f+fa;u=zeros(2*numnod,1);for i=1:numnod u2(1,i)= u(2*i-1);u2(2,i)= u(2*i);end nx1=2;ny1=10;I = Ly^3/12;

=

=

for i=1:(ny1+1)xjm(1,i)= Lx/2;xjm(2,i)=-(Ly)/ny1*(i-1)+Ly;yjm(i)=-(Ly/ny1)*(i-1)+Ly/2;stress11ex(i)=-q*(Lx-xjm(1,i))* yjm(i)/I;stress12ex(i)= q/(2*I)*(Ly^2/4-yjm(i)^2);end ind = 0;enorm=0;for gg=xjm ind = ind+1;gpos = gg(1:2);v = domain(gpos,x,dm,numnod);L = length(v);en = zeros(1,2*L);[phi,dphix,dphiy] = shape(gpos,dmax,x,v,dm);Bmat=zeros(3,2*L);for j=1:L Bmat(1:3,(2*j-1):2*j)= [dphix(j)0;0 dphiy(j);dphiy(j)dphix(j)];end for i=1:L en(2*i-1)= 2*v(i)-1;en(2*i)= 2*v(i);end

stress(1:3,ind)= Dmat*Bmat*u(en);%stressex(1,ind)=;% stressex(2,ind)= 0;

% stressex(3,ind)= 0;% err = stress(1:3,ind)-stressex(1:3,ind);% err2 = weight*jac*(0.5*(inv(Dmat)*err)'*(err));% enorm = enorm + err2;end %uex=zeros(2,numnod);I = Ly^3/12;ind4 = 0;for i=1:numnod if(x(2,i)==Ly/2)ind4=ind4+1;uex2(ind4)

= q/(6*young*I)*(3*nu*(x(2,i)-Ly/2)^2*(Lx-x(1,i))+(4+5*nu)*(Ly/2)^2*x(1,i)+(3*Lx-x(1,i))*x(1,i)^2);end figure hold on plot(x(1,(ny+1)/2:(ny+1):numnod),u2(2,(ny+1)/2:(ny+1):numnod),'r.');plot(x(1,(ny+1)/2:(ny+1):numnod),uex2,'-');%plot(xz,u2jy,'o');xlabel('x/m','fontweight','bold');ylabel('ux/m','fontweight','bold');legend('Uynode','Exact Solution');hold off % figure % hold on % plot(xjm(2,1:(ny1+1)),stress(1,1:ind),'r*');% plot(xjm(2,1:(ny1+1)),stress11ex(1,1:ind),'.-');% legend('EFG Solution','exact solution');

% % xlabel('y/m','fontweight','bold');% ylabel('Stress ','fontweight','bold');% % hold off % % figure % hold on % plot(xjm(2,1:(ny1+1)),stress(3,1:ind),'r*');% plot(xjm(2,1:(ny1+1)),stress12ex(1,1:ind),'.-');% legend('EFG Solution','exact solution');% xlabel('y/m','fontweight','bold');% ylabel('Stress ','fontweight','bold');hold off Toc 矩形区域内均匀节点布置：

解析解与无网格近似的比较：

第五篇：C++课程设计高斯消元法求线性代数方程组的解

河北工业大学计算机软件技术基础（VC）课程设计报告

学院 管理 班级 管理104班 姓名 杨立宝 __ 学号 101707____ 成绩 __ ____

一、题目：

求线性代数方程组的解（高斯消去法）（C13）

二、设计思路

1、总体设计

1）分析程序的功能

第一：编写输入程序，通过键盘先输入对应的已知量及函数的大小n和系数a[i]和得数b[i]。

第二：编写中间程序，通过函数的调用先定义线性代数方程，然后通过程序求出方程的梯形矩阵系数，并最终得出结果。

第三编写输出程序，输出最终结果。

2）系统总体结构：设计程序的组成模块，简述各模块功能。模块一：各函数的具体内容

A：三个输入函数，分别输入n，一维数组，二维数组。即输入已知量。

B：中间运算函数，计算是使得方程系数所成的矩阵成梯形矩阵，未知数的结果。即计算中间变量及结果。

C：最后输出函数，输出最后计算结果。模块二：各函数原型的声明 a写头文件。

b变量声明:存放输入数据的数组的声明，存放中间变量的数组的声明，存放运算结果的数组的声明。分别存放对应数据。c输入有关操作的文字

d函数调用，在运算中自动调用对应的函数解决对应问题。模块三:主函数

2、各功能模块的设计：说明各功能模块的实现方法 模块一：各个函数的声明，直接声明。

模块二：各函数都通过for循环来实现各个数组之间的基本运算。

3、设计中的主要困难及解决方案

在这部分论述设计中遇到的主要困难及解决方案。1）困难1 函数调用是怎么用？ 解决方案：

仔细阅读课本，以及同学之间的讨论，和老师的帮助。

4、你所设计的程序最终完成的功能 1）说明你编制的程序能完成的功能

输入线性代数的系数后，运行程序即可得到梯形矩阵和结果。2）准备的测试数据及运行结果

三、程序清单

如果是使用一个文件完成的程序，只需列出程序代码。如果是使用多文件完成的程序，首先说明程序中的代码存放在哪些文件中，说明文件名（例如：本程序包含first.cpp、second.cpp、third.cpp和all.h四个文件）；然后依次给出每个文件名及该文件清单，例如：

#include const N= 10;//设定矩阵大小范围 /* * 使用已经求出的x，向前计算x（供getx()调用）* double a[][] 系数矩阵 * double x[] 方程组解 * int i 解的序号 * int n 矩阵大小 * return 公式中需要的和 */ double getm(double a[N][N], double x[N], int i, int n){ double m = 0;int r;for(r=i+1;r

result = double(b[n-1]/a[n-1][n-1]);else //计算其他x值（对于公式中的求和部分，需要调用getm()函数）result = double((b[i]-getm(a,x,i,n))/a[i][i]);

return result;} void main(){ //double a[N][N] = {{2},{1,3,2},{1,2,2}};//double b[N] = {4,6,5};double a[N][N];//系数矩阵 double b[N];//右端项 double x[N];//方程组解 int i,j,k;int n=N;//矩阵大小 /*用户手工输入矩阵*/ cout<<“请输入系数矩阵的大小：”;cin>>n;cout<<“请连续输入矩阵值：”;for(i=0;i>a[i][j];} cout<<“请输入右端项：”;for(i=0;i>b[i];} /*显示原始矩阵*/ cout<<“n原始矩阵n”;for(i=0;i

/*进行高斯消去*/ for(j=0;j

/*显示处理后矩阵*/ cout<<“高斯消去后矩阵n”;for(i=0;i

/*回代方式解方程组*/ for(i=n-1;i>=0;i--){ x[i] = getx(a,b,x,i,n);}

/*显示方程组解*/ cout<<“nn方程组解n”;for(i=0;i

四、对该设计题目有何更完善的方案

1、对自己完成程序进行自我评价。设计过程中遇到很多问题，但经过和同学讨论，以及老师的解答和查阅资料加上我的努力最终写出了程序。

五、收获及心得体会

1、通过本次课程设计，自己在哪些方面的能力有所提高。

通过对该程序的编写，使我对数组，for循环以及函数的调用有了深刻的认识，锻炼了自己对c++更深一步的了解。

2、收获和心得体会。

明确的知道了数组，for循环以及函数的调用有的运用。明白了编程是一种严谨的工作，锻炼了自己的思维能力，并将努力培养自己严谨的思维。

日期： 2011年

6月

23日