當前位置:首頁 » 編程語言 » 三對角追趕法c語言
擴展閱讀
webinf下怎麼引入js 2023-08-31 21:54:13
堡壘機怎麼打開web 2023-08-31 21:54:11

三對角追趕法c語言

發布時間: 2022-10-28 14:42:59

A. 追趕法解三對角矩陣的問題

追趕法本質上就是不選主元的Gauss消去法, 只要所有順序主子式都非零就可以用
對角占優性可以保證順序主子式都非零, 同時保證數值穩定性
如果沒有對角占優性有可能會出現演算法中斷或者即使不中斷但精度也比較低的情況

B. 誰能幫我下載個csdn文檔, 實現追趕法求解三對角矩陣方程組的C++源代碼,這是題目

實現追趕法求解三對角矩陣方程組的C++源代碼.doc

網盤下載:

C. 求c語言編程用追趕法求解三對角方程組2x1+x2=1,x1+3x2+x3=2,x2+x3+x4=2,2x3+x4=0懸10分

//2x1+x2=1,x1+3x2+x3=2,x2+x3+x4=2,2x3+x4=0
#include "stdio.h"
#include "math.h"
#define MAXNUM 10 //變數數量
int array[MAXNUM][MAXNUM]={{2,1,0,0,1},{1,3,1,0,2},{0,1,1,1,2},{0,0,2,1,0}};

int unuse_result[MAXNUM];
int GaussFun(int equ,int var,int result[])
{
int i,j,k,col,num1,num2;
int max_r,ta,tb,gcdtemp,lcmtemp;
int temp,unuse_x_num,unuse_index;
col=0;
for (k=0;k<equ && col<var;k++,col++)
{
max_r=k;
for (i=k+1;i<equ;i++)
{
if (abs(array[i][col])>abs(array[max_r][col]))
{
max_r=i;
}
}
if (max_r!=k)
{
for (j=k;j<var+1;j++)
{
temp=array[k][j];
array[k][j]=array[max_r][j];
array[max_r][j]=temp;
}
}
if (array[k][col]==0)
{
k--;
continue;
}
for (i=k+1;i<equ;i++)
{
if (array[i][col]!=0)
{
num1=abs(array[i][col]);
num2=abs(array[k][col]);
while (num2!=0)
{
temp=num2;
num2=num1%num2;
num1=temp;
}
gcdtemp=num1;
lcmtemp=(abs(array[i][col])*abs(array[k][col]))/gcdtemp;
ta=lcmtemp/abs(array[i][col]);
tb=lcmtemp/abs(array[k][col]);
if (array[i][col]*array[k][col]<0)
{
tb=-tb;
}
for (j=col;j<var+1;j++)
{
array[i][j]=array[i][j]*ta-array[k][j]*tb;
}
}
}
}
for (i=k;i<equ;i++)
{
if (array[i][col]!=0)
{
return -1;
}
}
if (k<var)
{
for (i=k-1;i>=0;i--)
{
unuse_x_num=0;
for (j=0;j<var;j++)
{
if (array[i][j]!=0 && unuse_result[j])
{
unuse_x_num++;
unuse_index=j;
}
}
if (unuse_x_num>1)
{
continue;
}
temp=array[i][var];
for (j=0;j<var;j++)
{
if (array[i][j]!=0 && j!=unuse_index)
{
temp=array[i][j]*result[j];
}
}
result[unuse_index]=temp/array[i][unuse_index];
unuse_result[unuse_index]=0;
}
return var-k;
}
for (i=var-1;i>=0;i--)
{
temp=array[i][var];
for (j=i+1;j<var;j++)
{
if (array[i][j]!=0)
{
temp-=array[i][j]*result[j];
}
}
if (temp%array[i][i]!=0)
{
return -2;
}
result[i]=temp/array[i][i];
}
return 0;
}
void main()
{
int i,type;
int equnum,varnum;
int result[MAXNUM];
equnum=3;
varnum=3;
type=GaussFun(equnum,varnum,result);
if (type==-1)
{
printf("該方程無解!\n");
}
else if(type==-2)
{
printf("該方程有浮點數解,無整數解!\n");
}
else if (type>0)
{
printf("該方程有無窮多個解!自由變數數量為%d\n",type);
for (i=0;i<varnum;i++)
{
if (unuse_result[i])
{
printf("x%d 是不確定的\n",i+1);
}
else
{
printf("x%d:%d\n",i+1,result[i]);
}
}
}
else
{
printf("該方程的解為:\n");
for (i=0;i<varnum;i++)
{
printf("x%d=%d\n",i+1,result[i]);
}
}
}

D. 帶狀矩陣分解及追趕法

如果矩陣的非零元素分布在對角線元素兩側,例如,對角線下側有m1條次對角線為非零元素,對角線上側有m2條次對角為非零元素,矩陣的其他次對角線上元素均為零元素,這樣的矩陣,稱為帶寬為m1+m2+1的帶狀矩陣。連續性問題離散化後得到的線性方程組,其系數矩陣往往是帶狀矩陣。

若帶狀矩陣的各階主子式不為零,帶寬為m1+m2+1的矩陣A,可分解成帶寬為m1+1的下三角陣L和帶寬為m2+1的上三角陣U的乘積。

當m1=m2=1時,該帶狀矩陣稱為三對角矩陣。三對角矩陣是最簡單,然而又是常遇的帶狀矩陣。例如,在建立三次樣條插值函數以及求解微分方程數值解等問題中,都會遇到求解這樣的線性方程組Ax=d,其中系數矩陣A是一個階數較高的三對角方陣:

地球物理數據處理基礎

簡記為A=[ai,bi,cin1(i=1,2,…,n),這里n指n階方陣,1指半帶寬(簡稱帶寬)。右端項d=(d1,d2…,dnT。此時方程組Ax=d稱為三對角方程組,通過對三對角矩陣的直接分解來求解的方法稱為追趕法。

假設系數矩陣A=[ai,bi,cin1存在克洛脫(Crout)分解A=LU。由於矩陣A是帶寬為1的三對角矩陣,所以L,U均為二對角陣,即

地球物理數據處理基礎

利用矩陣乘法,比較A=LU兩端對應的元素,便可得如下確定矩陣L與U元素mi,li,ui的計算公式:

地球物理數據處理基礎

對三對角方程組的系數矩陣A進行分解以後,求解Ax=d的問題就轉化為求解兩個特殊的二對角方程組:

地球物理數據處理基礎

容易得到遞推的計算式分別為

地球物理數據處理基礎

人們常把下標由小到大求解 的過程形象化地稱為「追」的過程,而將下標由大到小求解 的過程稱為「趕」的過程。因此,上述方法通常稱為追趕法。

從上面的討論,我們看到追趕法的實質就是克洛脫分解法在求解三對角方程組中的應用。這時,由於矩陣A的特殊性,在求解運算時將系數矩陣中大量零元素撇開,僅在三條對角線上進行,從而簡化計算公式。

三對角線方程組的系數矩陣A=[ai,bi,cin1應滿足什麼條件才能進行克洛脫分解?下面給出一個充分條件:

在三對角矩陣A=[ai,bi,cin1中,若

(1)ai≠0(i=2,3,…n)且ci≠0(i=1,2,…,n-1),

(2)按行對角占優,即

則Α非奇異且存在唯一的Α=LU分解。

從上面的計算公式得知,用追趕法求解三對角方程組的計算量小,僅需5n-4次乘除法運算,且計算穩定、存儲量小。

[例]用追趕法求解三對角方程組

解:(1)先分解A=LU,求L和U的元素,

m2=a2=-1,m3=a3=-2,m4=a4=-3,l1=b1=2

l2=b2-m2u1=5/2,l3=b3-m3u2=12/5,l4=b4-m4u3=5/4

u1=c1/l1=-1/2,u2=c2/l2=-4/5,u3=c3/l3=-5/4

(2)求解Ly=d,即追的過程,

地球物理數據處理基礎

(3)求解Ux=y,即趕的過程,

x4=y4=2,x3=y3-u3x4=3,x2=y2-u2x3=4,x1=y1-u1x2=5

所以方程組的解為x=(5,4,3,2)T

E. 求c語言編程

// 程序7.5 追趕法求解三對角方程組

#include <stdio.h>
#include <conio.h>
#include <math.h>

#define MAX_n 100
#define PRECISION 0.0000001

//解輸出
void SulutionOutput(float x[],int n)
{
int i;
for(i=1;i<=n;++i)
printf("\nx[%d]=%f",i,x[i]);
}

//三對角方程組元素輸入
void TriDiagonalMatrixInput(float a[],float b[],float c[],float f[],int n)
{
int i;
printf("Input b[1],c[1],f[1]:");
scanf("%f %f %f",&b[1],&c[1],&f[1]);
for(i=2;i<n;++i)
{
printf("Input a[%d],b[%d],c[%d],f[%d]:",i,i,i,i);
scanf("%f %f %f %f",&a[i],&b[i],&c[i],&f[i]);
}
printf("Input a[%d],b[%d],f[%d]:",n,n,n);
scanf("%f %f %f",&a[n],&b[n],&f[n]);
}

//三對角方程組求解 — 追趕法
int Z_G_method(float a[],float b[],float c[],float f[],int n)
{
int i;
//
c[1]/=b[1];
for(i=2;i<n;++i)
c[i]/=(b[i]-a[i]*c[i-1]);
//
f[1]/=b[1];
for(i=2;i<=n;++i)
f[i]=(f[i]-a[i]*f[i-1])/(b[i]-a[i]*c[i-1]);
//
for(i=n-1;i>0;--i)
f[i]-=c[i]*f[i+1];
}

void main()
{
int n;
float a[MAX_n],b[MAX_n],c[MAX_n],f[MAX_n];
printf("\nInput n=");
scanf("%d",&n);
TriDiagonalMatrixInput(a,b,c,f,n);
Z_G_method(a,b,c,f,n);
SulutionOutput(f,n);
}

/*
運行實例:

Input n=4
Input b[1],c[1],f[1]:2 1 1
Input a[2],b[2],c[2],f[2]:1 3 1 2
Input a[3],b[3],c[3],f[3]:1 4 1 3
Input a[4],b[4],f[4]:1 5 4

x[1]=0.294118
x[2]=0.411765
x[3]=0.470588
x[4]=0.705882
*/

F. 用matalb編寫程序,用追趕法求解三對角線性方程組:

function x=Trid(a,b,c,d)
% 追趕法求解三對角的線性方程組 Ax=d
% b為主對角線元素,a,c分別為次對角線元素,d為右端項
% A=[ b1 c1
% a2 b2 c2
% ......
% a_(n-1) b_(n-1) c_(n-1)
% a_(n) b_(n) ]
% b=[b1...b_(n)]
% a=[0 a2...a_(n)]
% c=[c1...c_(n-1)]

n=length(b);

u(1)=b(1);
for i=2:n
l(i)=a(i)/u(i-1);
u(i)=b(i)-l(i)*c(i-1);
end

y(1)=d(1);
for i=2:n
y(i)=d(i)-l(i)*y(i-1);
end

x(n)=y(n)/u(n);
for i=n-1:-1:1
x(i)=(y(i)-c(i)*x(i+1))/u(i);
end

G. 我在visual c++ 2005中編的追趕法在debug下與release版本下為何運行結果不同

release經過的代碼優化.
你可以在項目編譯的配置里關掉一些的優化選項試試看.

H. 什麼是追趕法,其通用程序是什麼

function x=zhuiganfa
%首先說明:追趕法是適用於三對角矩陣的線性方程組求解的方法,並不適用於其他類型矩陣。
%定義三對角矩陣A的各組成單元。方程為Ax=d
% b為A的對角線元素(1~n),a為-1對角線元素(2~n),c為+1對角線元素(1~n-1)。
% A=[2 -1 0 0
% -1 3 -2 0
% 0 -2 4 -3
% 0 0 -3 5]
a=[0 -1 -2 -3];c=[-1 -2 -3];b=[2 3 4 5];d=[6 1 -2 1];
n=length(b);
u0=0;y0=0;a(1)=0;
%「追」的過程
L(1)=b(1)-a(1)*u0;
y(1)=(d(1)-y0*a(1))/L(1);
u(1)=c(1)/L(1);
for i=2:(n-1)
L(i)=b(i)-a(i)*u(i-1);
y(i)=(d(i)-y(i-1)*a(i))/L(i);
u(i)=c(i)/L(i);
end
L(n)=b(n)-a(n)*u(n-1);
y(n)=(d(n)-y(n-1)*a(n))/L(n);
%「趕」的過程
x(n)=y(n);
for i=(n-1):-1:1
x(i)=y(i)-u(i)*x(i+1);