數(shù)值分析-計算方法上機題集

1、實驗課程： 數(shù)值分析 專 業(yè)： 信息與計算科學 班 級： 10070241 學 號： 18 姓 名： 付強 中北大學理學院實驗一 函數(shù)差值方法【實驗目的】運用數(shù)值分析方法及相應的數(shù)學軟件、函數(shù)差值方法解決函數(shù)的近似求解。【實驗內(nèi)容】對于給定的一元函數(shù)的n+1個節(jié)點值（j=0，1,n）。試用拉格朗日公式求其差值多項式或分段二次拉格朗日插值多項式。數(shù)據(jù)如下：x0.40.550.650.800.951.05y0.410750.578150.696750.901.001.25382球五次拉格朗日插值多項式或分段三次插值多形式，并計算的值。（提示：結(jié)果為）【實驗所使用的儀器設(shè)備與軟件平臺】 電腦、VC+

2、6.0.【實驗方法與步驟】 首先，在理論層面上了結(jié)各種插值法的優(yōu)缺點，理解各種方法的插值原理。 實驗一，采取的是拉格朗日插值法，基本公式如下：其中，因此，實現(xiàn)此公式的關(guān)鍵步驟在于用C語言循環(huán)語句，實現(xiàn)此差值公式，然后運用數(shù)組輸入具體實驗數(shù)字，求的所需要的結(jié)果。實驗C語言代碼如下：#include <math.h> #include <stdio.h> #include <stdlib.h> #include<malloc.h> float lag(int n,float *x,float *y,float t) int i,j; float p,

3、s,la; p=1;la=0; for(i=0;i<n;i+) for(j=0;j<n;j+) if(j!=i) p=p*(t-*(x+j)/(*(x+i)-*(x+j); s=*(y+i)*p; p=1; la=la+s; return(la); void main() int m,n,k,i; float *x,*y,la,t; printf("請輸入給定點的數(shù)量n"); scanf("%d",&n); x=(float*)calloc(n,sizeof(float); if(x=NULL) printf("內(nèi)存分配失敗

4、"); exit(1); y=(float*)calloc(n,sizeof(float); if(y=NULL) printf("內(nèi)存分配失敗"); exit(1); x0=0.4;x1=0.55;x2=0.65;x3=0.80;x4=0.95;x5=1.05; y0=0.41075; y1=0.57815; y2=0.69675; y3=0.90; y4=1.00;y5=1.25382; for(i=0;i<n;i+) printf("x%d=%f ",i,*(x+i); printf("y%d=%fn",i,*(

5、y+i); for(i=0;i+) printf("請輸入t的值(退出則輸入0)"); scanf("%f",&t); if(t=0) printf("確認退出請輸入0，計算f(0)請按輸入任意數(shù)字值"); scanf("%d",&m); if(m=0)break; printf("請輸入需要使用的點的個數(shù)"); scanf("%d",&k); la=lag(k,x,y,t); printf("插值結(jié)果是:n"); printf(&

6、quot;f(%f)=%fn",t,la); free(x); free(y); 【實驗結(jié)果】【結(jié)果分析與討論】 實驗結(jié)果與實驗內(nèi)容中的提示結(jié)果完全一樣，甚至比實驗提示中所給出的精度都要高。實驗二 正交多項式與曲線擬合【實驗目的】用最小二乘法進行曲線擬合,并比較曲線擬合效果，探索擬合曲線的選擇與擬合精度間的關(guān)系?！緦嶒瀮?nèi)容】 1、 編寫出legendre、Chebyshev多項式的程序；2、 從隨機的數(shù)據(jù)中找出其規(guī)律性，給出其近似表達式的問題，在生產(chǎn)實踐和科學實驗中大量存在，通常利用數(shù)據(jù)的最小二乘法求得擬合曲線。在某冶煉過程中，根據(jù)統(tǒng)計數(shù)據(jù)的含碳量與時間關(guān)系，試求含碳量與時間的擬合曲

7、線。051015202530354045505501.272.162.863.443.874.154.374.514.584.024.6411111111111【實驗所使用的儀器設(shè)備與軟件平臺】 電腦、VC+6.0.、matlab【實驗方法與步驟】一、按照題目要求，首先對數(shù)據(jù)進行一次擬合，代碼如下：#include <math.h> #include <stdio.h> #include <stdlib.h> #include<malloc.h> float average(int n,float *x) int i; float av; av=

8、0; for(i=0;i<n;i+) av+=*(x+i); av=av/n; return(av); /平方和 float spfh(int n,float *x) int i; float a,b; a=0; for(i=0;i<n;i+) a+=(*(x+i)*(*(x+i); return(a); /和平方 float shpf(int n,float *x) int i; float a,b; a=0; for(i=0;i<n;i+) a=a+*(x+i); b=a*a/n; return(b); /兩數(shù)先相乘，再相加 float dcj(int n,float *

9、x,float *y) int i; float a; a=0; for(i=0;i<n;i+) a+=(*(x+i)*(*(y+i); return(a); /兩數(shù)先相加，再相乘 float djc(int n,float *x,float *y) int i; float a=0,b=0; for(i=0;i<n;i+) a=a+*(x+i); b=b+*(y+i); a=a*b/n; return(a); /系數(shù)a float xsa(int n,float *x,float *y) float a,b,c,d,e; a=spfh(n,x); b=shpf(n,x); c=d

10、cj(n,x,y); d=djc(n,x,y); e=(c-d)/(a-b); /printf("%f %f %f %f",a,b,c,d); return(e); float he(int n,float *y) int i; float a; a=0; for(i=0;i<n;i+) a=a+*(y+i); return(a); float xsb(int n,float *x,float *y,float a) float b,c,d; b=he(n,y); c=he(n,x); d=(b-a*c)/n; return(d); int main() int n,

11、i;float *x,*y,a,b; printf("請輸入將要輸入的有效數(shù)值組數(shù)n的值:"); scanf("%d",&n); x=(float*)calloc(n,sizeof(float); if(x=NULL) printf("內(nèi)存分配失敗"); exit(1); y=(float*)calloc(n,sizeof(float); if(y=NULL) printf("內(nèi)存分配失敗"); exit(1); printf("請輸入x的值n"); for(i=0;i<n;i+)

12、scanf("%f",x+i); printf("請輸入y的值,請注意與x的值一一對應:n"); for(i=0;i<n;i+)scanf("%f",y+i); for(i=0;i<n;i+) printf("x%d=%3.2f ",i,*(x+i); printf("y%d=%3.2fn",i,*(y+i); a=xsa(n,x,y); b=xsb(n,x,y,a); printf("經(jīng)最小二乘法擬合得到的一元線性方程為:n"); printf("f(

13、x)=%3.2fx+%3.2fn",a,b); 二、經(jīng)過matlab進行曲線擬合可知，該組數(shù)據(jù)近似的服從3次曲線，因此下面對所給數(shù)據(jù)進行3次最小二乘法擬合，代碼如下：#include <stdio.h>#include <conio.h>#include <math.h>#include <process.h>#define N 12/N個點#define T 3 /T次擬合#define W 1/權(quán)函數(shù)#define PRECISION 0.00001float pow_n(float a,int n)int i;if(n=0)ret

14、urn(1);float res=a;for(i=1;i<n;i+)res*=a;return(res);void mutiple(float aN,float bT+1,float cT+1)float res=0;int i,j,k;for(i=0;i<T+1;i+)for(j=0;j<T+1;j+)res=0;for(k=0;k<N;k+)res+=aik*bkj;cij=res;void matrix_trans(float aT+1,float bN)int i,j;for(i=0;i<N;i+)for(j=0;j<T+1;j+)bji=aij;v

15、oid init(float x_y2,int n)int i;printf("請輸入%d個已知點：n",N);for(i=0;i<n;i+)printf("(x%d y%d):",i,i);scanf("%f %f",&x_yi0,&x_yi1);void get_A(float matrix_AT+1,float x_y2,int n)int i,j;for(i=0;i<N;i+)for(j=0;j<T+1;j+)matrix_Aij=W*pow_n(x_yi0,j);void print_arr

16、ay(float arrayT+1,int n)int i,j;for(i=0;i<n;i+)for(j=0;j<T+1;j+)printf("%-g",arrayij);printf("n");void convert(float arguT+2,int n)int i,j,k,p,t;float rate,temp;for(i=1;i<n;i+)for(j=i;j<n;j+)if(argui-1i-1=0)for(p=i;p<n;p+)if(argupi-1!=0)break;if(p=n)printf("方程

17、組無解!n");exit(0);for(t=0;t<n+1;t+)temp=argui-1t;argui-1t=argupt;argupt=temp;rate=arguji-1/argui-1i-1;for(k=i-1;k<n+1;k+)argujk-=argui-1k*rate;if(fabs(argujk)<=PRECISION)argujk=0;void compute(float arguT+2,int n,float root)int i,j;float temp;for(i=n-1;i>=0;i-)temp=arguin;for(j=n-1;j&g

18、t;i;j-)temp-=arguij*rootj;rooti=temp/arguii;void get_y(float trans_AN,float x_y2,float y,int n)int i,j;float temp;for(i=0;i<n;i+)temp=0;for(j=0;j<N;j+)temp+=trans_Aij*x_yj1;yi=temp;void cons_formula(float coef_AT+1,float y,float coef_formT+2)int i,j;for(i=0;i<T+1;i+)for(j=0;j<T+2;j+)if(j

19、=T+1)coef_formij=yi;elsecoef_formij=coef_Aij;void print_root(float a,int n)int i,j;printf("%d個點的%d次擬合的多項式系數(shù)為:n",N,T);for(i=0;i<n;i+)printf("a%d=%g,",i+1,ai);printf("n");printf("擬合曲線方程為:ny(x)=%g",a0);for(i=1;i<n;i+)printf(" + %g",ai);for(j=0;j&l

20、t;i;j+)printf("*X");printf("n");void process()float x_yN2,matrix_ANT+1,trans_AT+1N,coef_AT+1T+1,coef_formuT+1T+2,yT+1,aT+1;init(x_y,N);get_A(matrix_A,x_y,N);printf("矩陣A為：n");print_array(matrix_A,N);matrix_trans(matrix_A,trans_A);mutiple(trans_A,matrix_A,coef_A);printf(&

21、quot;法矩陣為：n");print_array(coef_A,T+1);get_y(trans_A,x_y,y,T+1);cons_formula(coef_A,y,coef_formu);convert(coef_formu,T+1);compute(coef_formu,T+1,a);print_root(a,T+1);void main()process();【實驗結(jié)果】【結(jié)果分析與討論】對一次擬合做matlab檢驗，結(jié)果如下： 可從圖像得一次擬合是錯誤的，于是進行三次擬合，結(jié)果如下： 因此，得出上述數(shù)據(jù)應進行三次最小二乘法擬合才能得到比較精確的解。實驗三 數(shù)值積分與數(shù)值微

22、分【實驗目的】選用復合梯形公式，復合辛普森公式，算法，高斯算法計算所給出的積分的近似解?！緦嶒瀮?nèi)容】選用復合梯形公式，復合辛普森公式，算法，高斯算法計算（1）（2）（3）【實驗所使用的儀器設(shè)備與軟件平臺】電腦、VC+6.0?！緦嶒灧椒ㄅc步驟】一、運用算法，對（1）進行數(shù)值法幾分計算代碼如下：#include <stdio.h>#include <math.h>static double T2020; double f(double x)return pow(4-pow(sin(x),2),0.5); double fuhe(double a,double b,int n

23、)double sum=0,h=0,k=0;int i;if(n=0)h=b-a;return h*(f(a)+f(b)/2.0;elsek=pow(2,n);h=(b-a)/k;for(i=0;i<k;i+)sum+=(f(a+(b-a)*i/k)+f(a+(b-a)*(i+1)/k)*h/2.0;return sum; double romberg(double a,double b,double e)int j=0,k=0,i=2;T01=fuhe(a,b,0);T11=fuhe(a,b,1);T02=(4*T11-T01)/3.0;for(;fabs(T0i-T0i-1)>

24、=e;i+)Ti1=fuhe(a,b,i+1);j=2;for(k=i-1;k>=0;k-)Tkj=(pow(4,j-1)*Tk+1j-1-Tkj-1)/(pow(4,j-1)-1);j+;return T0i; int main()float a,b,e;double result;printf("請輸入積分下限a:n");scanf("%f",&a);printf("請輸入積分上限b:n");scanf("%f",&b);printf("輸入e的近似值:n");scan

25、f("%f",&e);result=romberg(a,b,e);printf("Based on your input parameters the results is:n%fnn",result);scanf("%f",&e);二、運用復合辛普森算法，對（1）進行數(shù)值法幾分計算代碼如下：#include<stdio.h>#include<math.h>double f(double x); double zdyf(double x); double ff(double a,double b

26、,double eps);double ff(double a,double b,double eps,double (*f)(double);void main()char str1;double a,b,eps,t; /a,b內(nèi),當前精度eps,接受結(jié)果到t;double (*fff)(double);a=0.0;b=0.25;eps=0.0000001;t=ff(a,b,eps);printf("被積函數(shù)：n");printf(" pow(4-pow(sin(x),2),0.5)n");printf("積分運算：n");prin

27、tf(" 區(qū)間%3.2f,%3.2fn",a,b);printf("積分結(jié)果:n");printf(" t=%en",t);printf("n");gets(str);/*-自定義被積函數(shù)-*/double f(double x)double y;y=pow(4-pow(sin(x),2),0.5);return(y);/*-積函數(shù)（辛普森）在區(qū)間a,b積f(x)-*/double ff(double a,double b,double eps) int n,k; double h,t1,t2,s1,s2,ep,p

28、,x; n=1; /當前子節(jié)點個數(shù)-1，不包括兩端a,b h=b-a; /當前區(qū)間長度 t1=h*(f(a)+f(b)/2.0; /計算原始梯形面積 s1=t1; /將原始梯形面積->s1 ep=eps+1.0; /誤差初始化/*-辛普森算法-*/while(ep>=eps) /循環(huán)比較(當前ep)和(標準eps) p=0.0; /累加暫存量初始化 for(k=0;k<=n-1;k+) /循環(huán)求和實現(xiàn) x=a+(k+0.5)*h; /節(jié)點定位 p+=f(x); /節(jié)點函數(shù)值累加 t2=(t1+h*p)/2.0; /用(上一個t1)構(gòu)造(當前t2); s2=(4.0*t2-t1

29、)/3.0; /用 (上一個t1)和(當前t2) 構(gòu)造(當前s2) ep=fabs(s2-s1); /用 (當前s2)和(上一個s1) 構(gòu)造(當前ep) t1=t2; /為構(gòu)造下一個t準備 s1=s2; /為構(gòu)造下一個s準備 n+=n; /節(jié)點個數(shù)倍增 h/=2.0; /區(qū)間二分化return(s2); /返回(當前s2)double ff(double a,double b,double eps,double (*f)(double) int n,k; double h,t1,t2,s1,s2,ep,p,x; n=1; /當前子節(jié)點個數(shù)-1，不包括兩端a,b h=b-a; /當前區(qū)間長度 t

30、1=h*(*f)(a)+(*f)(b)/2.0; /計算原始梯形面積 s1=t1; /將原始梯形面積->s1 ep=eps+1.0; /誤差初始化/*-辛普森算法-*/while(ep>=eps) /循環(huán)比較(當前ep)和(標準eps) p=0.0; /累加暫存量初始化 for(k=0;k<=n-1;k+) /循環(huán)求和實現(xiàn) x=a+(k+0.5)*h; /節(jié)點定位 p+=(*f)(x); /節(jié)點函數(shù)值累加 t2=(t1+h*p)/2.0; /用(上一個t1)構(gòu)造(當前t2); s2=(4.0*t2-t1)/3.0; /用 (上一個t1)和(當前t2) 構(gòu)造(當前s2) ep=

31、fabs(s2-s1); /用 (當前s2)和(上一個s1) 構(gòu)造(當前ep) t1=t2; /為構(gòu)造下一個t準備 s1=s2; /為構(gòu)造下一個s準備 n+=n; /節(jié)點個數(shù)倍增 h/=2.0; /區(qū)間二分化return(s2); /返回(當前s2)【實驗結(jié)果】算法 復合辛普森算法 【結(jié)果分析與討論】經(jīng)過對比兩種方法的實驗結(jié)果，可以知道兩種方法的計算結(jié)果正確，并且具有一定的代數(shù)精度。實驗四 常微分方程初值問題數(shù)值解法【實驗目的】利用Euler法，改進Euler法，Rung-Kutta方法求常微分方程數(shù)值解【實驗內(nèi)容】科學計算中經(jīng)常遇到微分方程（組）初值問題，需要利用Euler法，改進Euler

32、法，Rung-Kutta方法求其數(shù)值解，諸如以下問題：分別取h=0.1,0.2,0.4時數(shù)值解。初值問題的精確解。【實驗所使用的儀器設(shè)備與軟件平臺】電腦、VC+6.0?！緦嶒灧椒ㄅc步驟】一、利用Euler法求數(shù)值解，代碼如下：#include <iostream>using namespace std;#include <math.h>float f(float x, float y);void init();float h, /步長 x00, /x的取值范圍 xn, / x0< x < xn y00; /初值float x21;/結(jié)果float y21;v

33、oid Euler() init(); x0 = x00; y0 = y00; for (int i=1; i<=(xn-x00)/h+1; i+) xi = xi-1 + h; yi = yi-1 + h * f(xi-1,yi-1); float f(float x, float y) return 4*x/y-x*y; /return y-2*x/y;/初始化初值和x的范圍void init() x00=0; xn=2; y00=1;/求常微分方程的標準解float pf(float x) return sqrt(4+5*exp(-x*x);int main() h=0.1; Eu

34、ler(); cout << "步長 h = " << h << endl; for (int i=1; i<=(xn-x00)/h+1; i+) cout << "x" << i << "=" << xi << "t" << "y" << i << "=" << yi << "t" <&l

35、t; "y("<< xi << ")=" << pf(xi) << endl; h=0.2; Euler(); cout << "步長 h = " << h << endl; for (i=1; i<=(xn-x00)/h+1; i+) cout << "x" << i << "=" << xi << "t" <<

36、"y" << i << "=" << yi << "t" << "y("<< xi << ")=" << pf(xi) << endl; h=0.4; Euler(); cout << "步長 h = " << h << endl; for (i=1; i<=(xn-x00)/h+1; i+) cout << &qu

37、ot;x" << i << "=" << xi << "t" << "y" << i << "=" << yi << "t" << "y("<< xi << ")=" << pf(xi) << endl; return 0;二、用四階龍格庫塔法求數(shù)值解，matlab代碼如下（以h=

38、0.1為例）：function varargout=saxplaxliu(varargin)clc,clearx0=0;xn=1.0;y0=1;h=0.1;y,x=lgkt4j(x0,xn,y0,h); n=length(x);fprintf(' i x(i) y(i)n');for i=1:n fprintf('%2d %f %fn',i,x(i),y(i);endfunction z=f(x,y) z=4*x/y-x*y;function y,x=lgkt4j(x0,xn,y0,h) x=x0:h:xn;n=length(x);y1=x;y1(1)=y0;f

39、or i=1:n-1 K1=f(x(i),y1(i); K2=f(x(i)+h/2,y1(i)+h/2*K1); K3=f(x(i)+h/2,y1(i)+h/2*K2); K4=f(x(i)+h,y1(i)+h*K3); y1(i+1)=y1(i)+h/6*(K1+2*K2+2*K3+K4);end y=y1;【實驗結(jié)果】一、C語言Euler法求數(shù)值解運行結(jié)果：二、matlab四階龍格庫塔運行結(jié)果：【結(jié)果分析與討論】通過比較C語言、matlab兩種不同運行環(huán)境下，兩種不同的算法所得出的結(jié)果是一樣的。因此實驗方法得當，結(jié)論正確。對于歐拉法對于歐拉法，步長越小，精度越高，而產(chǎn)生的誤差越小?？傮w來說

40、，歐拉法的優(yōu)點是形式簡單，計算方便，缺點是總的運算精度比較低。而且隨著x的增大，誤差值也越來越大。根據(jù)歐拉公式的截斷誤差計算，歐拉法是一階方法。對于龍格庫塔法，優(yōu)點是精度更高，同樣的步長下精度比歐拉法高的多，誤差更小，截斷誤差為O(h5)，故龍格庫塔法是四階方法。缺點是每步都要計算四次微分值。實驗五 非線性方程求根【實驗目的】1、 編制一個程序進行運算，最后打印出每種迭代格式的斂散情況2、 用事后誤差估計來控制迭代次數(shù)，并且打印出迭代的次數(shù)3、 初始值的選取對迭代收斂有何影響4、 分析迭代收斂和發(fā)散的原因【實驗內(nèi)容】設(shè)方程有三個實根，現(xiàn)采用下面四種不同計算格式，求的根。1、（，則從滿足收斂條件

41、的角度來說，則只能求）2、（，則從滿足收斂條件的角度來說，則只能求）3、（，則從滿足收斂條件的角度來說，則只能求，）4、（，則從滿足收斂條件的角度來說，且，則只能求）【實驗所使用的儀器設(shè)備與軟件平臺】電腦、VC+6.0?！緦嶒灧椒ㄅc步驟】1、#include<stdio.h>#include<math.h>void main() double x0, x, e; int i, n, flag; flag = 0; e = 0.0001; n = 10000; x0 = -1; for(i=1;i<n;i+) x = (3 * x0 + 1) / x0 / x0;

42、if(fabs(x - x0) < e) printf("x = %lf, error < %lf, 迭代次數(shù) = %d",x,e,i); flag = 1; break; x0 = x; if (flag=0) printf("看來不收斂 x = %lf, error >= %lf, 迭代次數(shù) = %d", x,e,n); 2、#include<stdio.h>#include<math.h>void main() double x0, x, e; int i, n, flag; flag = 0; e = 0

43、.0001; n = 10000; x0 = 0; for(i=1;i<n;i+) x = (x0*x0*x0-1)/3.0; if(fabs(x - x0) < e) printf("x = %lf, error < %lf, 迭代次數(shù) = %d",x,e,i); flag = 1; break; x0 = x; if (flag=0) printf("看來不收斂 x = %lf, error >= %lf, 迭代次數(shù) = %d", x,e,n);3、#include<stdio.h>#include<math

44、.h>void main() double x0, x, e; int i, n, flag; flag = 0; e = 0.0001; n = 10000; x0 = 2; for(i=1;i<n;i+) x = pow(3*x0+1,0.3333333); if(fabs(x - x0) < e) printf("x = %lf, error < %lf, 迭代次數(shù) = %d",x,e,i); flag = 1; break; x0 = x; if (flag=0) printf("看來不收斂 x = %lf, error >=

45、 %lf, 迭代次數(shù) = %d", x,e,n); 4、#include<stdio.h>#include<math.h>void main() double x0, x, e; int i, n, flag; flag = 0; e = 0.0001; n = 10000; x0 = 2; for(i=1;i<n;i+) x = pow(3+pow(x,-1),0.5); if(fabs(x - x0) < e) printf("x = %lf, error < %lf, 迭代次數(shù) = %d",x,e,i); flag

46、= 1; break; x0 = x; if (flag=0) printf("看來不收斂 x = %lf, error >= %lf, 迭代次數(shù) = %d", x,e,n); 【實驗結(jié)果】1、2、3、4、【結(jié)果分析與討論】實驗所得結(jié)果與元命題給出的結(jié)果一致，唯一的缺點就是，在不收斂的點處求不出近似值，也運行不出結(jié)果。實驗六 線性方程組的直接解法【實驗目的】1、通過該課題的實驗，體會模塊化結(jié)構(gòu)程序設(shè)計方法的優(yōu)點；2、運用所學的計算方法，解決各類線性方程組的直接算法；3、通過三對角線性方程組的解法，體會稀疏線性方程組解法的優(yōu)點。【實驗內(nèi)容】給出下列幾個不同類型的線性方程

47、組1、（用高斯列主元消去法或直接LU分解法）2、 （用追趕法）【實驗所使用的儀器設(shè)備與軟件平臺】電腦、VC+6.0?！緦嶒灧椒ㄅc步驟】一、 高斯列主元消去法求1、2的解，代碼如下：#include<stdio.h>#include <math.h>#include<conio.h>#define N 20int main() int n,i,j,k; int mi,tmp,mx; float aNN,bN,xN; printf("/nInput n:"); scanf("%d",&n); if(n>N)

48、printf("The input n should in(0,N)!/n"); getch(); return 1; if(n<=0) printf("The input n should in(0,N)!/n"); getch(); return 1; printf("Now input a(i,j),i,j=0.%d:/n",n-1); for(i=0;i<n;i+) for(j=0;j<n;j+) scanf("%f",&aij); printf("Now input b(i),i,j=0.%d:/n",n-1); for(i=0;i<n;i+) scanf("%f",&bi); for(i=0;i<n-2;i+) for(j=i+1,mi=i,mx=fabs(aij);j<n-1;j+) if(fabs(aji)>mx) mi=j; mx=fabs(aji); if(i<mi) tmp=bi;bi=bmi;bmi=tmp; for(j=i;j<n;j+) tmp=aij; aij=amij; amij=tmp; for(j=i+1;j<n