实验一 直接法解线方程组.doc

实验一 直接法解线性方程组/1-1 Gauss列主元消去法【源代码】：#include #include #include#define eps 1e-9double *A,*X;int N;int FindMax(int k);void ExchangeRow(int k,int r);void breakup();void back();int main() int i,j; printf(请输入增广矩阵(A,B)的行数:t); scanf(%d,&N); A=(double*)calloc(N*(N+1),sizeof(double); X=(double*)calloc(N,sizeof(double); if(N256|N=0) printf(输入的数字不再范围之内!); printf(n); return 0; else /A=(double*)calloc(N*(N+1),sizeof(double); printf(请输入待求解方程组的增广矩阵(%d行 %d列):n,N,N+1); for(i=0;iN;i+) for(j=0;j(N+1);j+) scanf(%lf,&Aj+i*(N+1); printf(n增广矩阵(A,B)如下:n); for(i=0;iN;i+) for(j=0;j(N+1);j+) printf(%lft,Aj+i*(N+1); printf(n); breakup(); /上三角变换 back(); /回代函数 printf(n方程组的解为:n); for(i=0;iN;i+) printf(X(%d)=t%lfn,i+1,Xi); free(A); free(X); exit(0);/在列向量中寻找绝对值最大的项，并返回该项的标号int FindMax(int k) int i,r; double max=0.0; for(i=k;imax) max=fabs(Ai*(N+1)+k); r=i; /r保存最大项的行标号 return r;/换行void ExchangeRow(int k,int r) int j; double temp=0.0; for(j=k;jN+1;j+) temp=Ak*(N+1)+j; Ak*(N+1)+j=Ar*(N+1)+j; Ar*(N+1)+j=temp; /上三角变换，A为增广矩阵的指针，N为矩阵的行数。void breakup() int p,j,i,r; double m=0.0,det; for(i=0;iN-1;i+) r=FindMax(i); /找出该列最大项的标号 if(i!=r) ExchangeRow(i,r); /交换i行和r行 if(fabs(Ai*(N+1)+i)=0|fabs(Ai*(N+1)+i)=eps) printf(矩阵是一个奇异矩阵,没有唯一解!); break; /消去P元素以下的p列内容。 for(j=i+1;jN;j+) m=Aj*(N+1)+i/Ai*(N+1)+i; for(p=i;pN+1;p+) Aj*(N+1)+p=Aj*(N+1)+p-m*Ai*(N+1)+p; printf(n增广矩阵高斯列主元消去后的矩阵如下:n); for(i=0;iN;i+) for(j=0;j(N+1);j+) printf(%lft,Aj+i*(N+1); printf(n); det=1; for(i=0;i=0;k-) total=0.0; for(i=k+1;iN;i+) total=total+Ak*(N+1)+i*Xi; Xk=(Ak*(N+1)+N-total)/Ak*(N+1)+k; 【运行结果】：/1-2 追赶法【源代码】：#include #include #include #define N 1000/#define eps 1e-9int n;double ANN=0,fN,XN;double aN,bN,cN,yN;void divide() /Crout分解 int i; a0=A00; for(i=0;in;i+) ci=Ai+1i; /a0=A00; bi=Aii+1/ai; ai+1=Ai+1i+1-ci*bi; void back() /回代求解 int i; y0=f0/a0; for(i=1;i=0;i-) Xi=yi-bi*Xi+1;int main() int i,j; printf(输入矩阵A的阶数:); scanf(%d,&n);printf(输入矩阵A%d%d:n,n,n); scanf(%lf %lf,&A00,&A01); for(i=1;in-1;i+) scanf(%lf %lf %lf,&Aii-1,&Aii,&Aii+1); scanf(%lf %lf,&Aii-1,&Aii); printf(n输入右端项f: ); for(i=0;in;i+) scanf(%lf,&fi); divide(); back(); for(i=0;in;i+) printf(X%d=%lfn,i,Xi); return 0;【运行结果】：实验二 插值方法/2-1 Lagrange插值法【源程序】：#include #include #include double lagrange(double *x,double *y,double xx,int n) int i,j; double *a,yy=0; a=(double*)malloc(n*sizeof(double); for(i=0;in;i+) ai=yi; for(j=0;jn;j+) if(j!=i) ai*=(xx-xj)/(xi-xj); yy+=ai; free(a); return yy;int main() int i,n; double x100,y100,xx,yy; printf(Input n:); scanf(%d,&n); printf(n数组x:); for(i=0;in;i+) scanf(%lf,&xi); printf(n数组y:); for(i=0;in;i+) scanf(%lf,&yi); printf(nInpute xx:); scanf(%lf,&xx); yy=lagrange(x,y,xx,n); printf(x=%lfty=%lfn,xx,yy); return 0;/2-2 Newton插值法（画图）【源程序】：import java.awt.*; import java.awt.event.*;import java.awt.geom.*;import javax.swing.*;public class Lagrange extends JFrame JLabel ln=new JLabel(节点数n:), lxx=new JLabel(待求x:), lf=new JLabel(f(x):);JTextField tn=new JTextField(5), txx=new JTextField(20), tf=new JTextField(20);JButton b1=new JButton(确定), bpaint=new JButton(画图), b2=new JButton(求值);JTextField textx=new JTextField100, texty=new JTextField100;int i=0,n=0, width=100, height=100;public Lagrange()setTitle(插值方法-画图); setLayout(null);setBounds(100, 100, 500, 375);ln.setBounds(30, 30, 80, 25); /节点数ntn.setBounds(100, 30, 40, 25);b1.setBounds(160, 30, 60, 25); /按钮确定lxx.setBounds(30, 150, 80, 25); /待求xtxx.setBounds(100, 150, 100, 25);b2.setBounds(210, 150, 60, 25);/按钮求值lf.setBounds(30, 200, 80, 15);/f(x)tf.setBounds(100, 200, 100, 25);bpaint.setBounds(210, 200, 60, 25);/按钮画图ln.setFont(new Font(, 1, 15);tn.setFont(new Font(, 1, 15);b1.setFont(new Font(,1, 12);bpaint.setFont(new Font(, 1, 12);lxx.setFont(new Font(, 1, 15);txx.setFont(new Font(, 1, 15);b2.setFont(new Font(,1, 12);lf.setFont(new Font(,1, 15);tf.setFont(new Font(, 1, 15);add(ln); add(tn); add(b1); add(bpaint);add(lxx); add(txx); add(b2);add(lf); add(tf);add(bpaint);for(i=0;i15;i+)textxi=new JTextField();textxi.setBounds(30+i*80, 80, 80, 25);textxi.setFont(new Font(,1,12);add(textxi);textxi.setVisible(false);for(i=0;i15;i+)textyi=new JTextField();textyi.setBounds(30+i*80, 105, 80, 25);textyi.setFont(new Font(,1,12);add(textyi);textyi.setVisible(false);b1.addActionListener(new ActionListener()public void actionPerformed(ActionEvent e)if(!tn.getText().equals()n=Integer.parseInt(tn.getText();for(i=0;in;i+)textxi.setVisible(false);textyi.setVisible(false);for(i=0;in;i+)textxi.setText();textxi.setVisible(true);for(i=0;in;i+)textyi.setText();textyi.setVisible(true););bpaint.addActionListener(new ActionListener()public void actionPerformed(ActionEvent e)paint(););b2.addActionListener(new ActionListener()public void actionPerformed(ActionEvent e)Double res;res = function(Double.parseDouble(txx.getText();tf.setText(res.toString(););setVisible(true);double x=new double100, y=new double100, l=new double100;public double function(double xx) /计算插入的 xx的函数值 f(x) int j=0;Double yy=0.0; for(i=0;in;i+)if(!textxi.getText().equals()xi=Double.parseDouble(textxi.getText();if(!textyi.getText().equals()yi=Double.parseDouble(textyi.getText();for(i=0; in; i+)li=1; for(j=0; jn; j+) if(j!=i) li*=(xx-xj)/(xi-xj); yy+=yi*li;return yy;public void paint() /画函数图线class painter extends JFramepublic painter()super(绘制函数曲线);setLayout(null);setBackground(Color.PINK);this.setSize(450,450);setLocation(400, 200);setVisible(true);Label l=new Label(0,0);/坐标原点l.setFont(new Font(,15,15);l.setBounds(195, 180, 35, 20);add(l);public void paint(Graphics g)Line2D.Double line;Graphics2D g2 = (Graphics2D) g;g2.setFont(new Font(,30,30);g.drawLine(20, 200, 400, 200); /画x轴 g.drawLine(390, 190, 400, 200);/箭头g.drawLine(390, 210, 400, 200);g.drawLine(200, 30, 200, 400);/画y轴 g.drawLine(190, 40, 200, 30);/箭头g.drawLine(210, 40, 200, 30);for(double i=-200;i=200;i+=1e-3)double y=(double)function(i);line = new Line2D.Double(i*10+200, 200-y*10, i*10+200, 200-y*10);g2.draw(line);painter p=new painter();Graphics g=p.getGraphics();g.setFont(new Font(,30,30);p.paint(g);public static void main(String args)new Lagrange();【运行结果】： 实验三 数值计算/3-1 变步长梯形法求定积分【源程序】：#include #include using namespace std;double f(double x)/预先输入的待积分函数 double y; if(x=0) y=1; else y=sin(x)/x; return(y);double fun(double a,double b,double eps)/变步长梯形法 int n,i; double h,sum,t1000; n=1; t1=(b-a)*(f(a)+f(b)/2; do h=(b-a)/n; sum=0; for(i=0;ieps); return(tn);int main() double a,b,eps,result; cout abeps; coutendl; result=fun(a,b,eps); cout积分结果：resultendl; return 0;【运行结果】：/3-2 Romberge求定积分【源程序】：#include#includeusing namespace std;double T10001000;double fun(const double x);/求被积函数并返回double Romberge(const double a,const double b,const double eps); /求定积分函数(a,b,eps)int main()double a,b,eps,result;coutabeps;coutendl;result=Romberge(a,b,eps);cout所求积分：resultendl;return 0;double fun(const double x)if(x=0) return 1;return(sin(x)/x);double Romberge(const double a,const double b,const double eps)int l=0;T00=(b-a)*(fun(a)+fun(b)/2;dol+;int i,m,k,n=pow(2,l-1); /n为区间等分数double h=(b-a)/n; /h为步长 double temp=0; for(i=0;in;i+)temp+=fun(a+(i+0.5)*h);T0l=T0l-1/2+h/2*temp;for(m=1;m=l;+m) int t=l-m+1; for(k=1;keps);return (Tl0);/返回符合精度的值【运行结果】：实验四 常微分方程的数值解/4-1 改进欧拉法【源程序】：#include #include double Y;double f(double x,double y) return (-x*y*y);int main() double a,b,h,y0,temp,xx; int i; printf(a,b,h,y0: ); scanf(%lf%lf%lf%lf,&a,&b,&h,&y0); Y=y0; for(i=0,xx=a;i=(b-a)/h;i+,xx+=h) temp=Y+h*f(xx,Y); Y=Y+h*(f(xx,Y)+f(xx+h,temp)/2; printf(x=%lfty%d=%lfn,xx+h,i,Y); return 0;【运行结果】：/4-2 四阶R-K法【源程序】：#include#include#include#include double *X,*Y;int N;double function(double x,double y) return (-x*y*y); int main() double y0,a,b,h; int i; printf(Input a,b,h,y0: ); scanf(%lf%lf%lf%lf,&a,&b,&h,&y0); N=(b-a)/h; X=(double*)malloc(N+1)*sizeof(double); Y=(double*)malloc(N+1)*sizeof(double); double k1,k2,k3,k4;X0=a; Y0=y0;for(i=0;iN;i+) Xi+1=Xi+h; k1=function(Xi,Yi); k2=function(Xi+h/2,Yi+h*k1/2); k3=function(Xi+h/2,Yi+h*k2/2); k4=function(Xi+h,Yi+h*k3); Yi+1=Yi+h*(k1+2*k2+2*k3+k4)/6; for(i=0;iN+1;i+) printf(x%d=%lfty%d=%lfn,i,Xi,i,Yi); return 0;【运行结果】：实验五 迭代法解线性(非线性)方程/5-1 G-S迭代法解线性方程组【源程序】：#include #include#include#includeusing namespace std;const int N=100;double ANN,bN,X0N,X1N;int main() int i,j,n; double eps; cout neps; coutendl; printf(输入系数矩阵A%d%d:n,n,n); for(i=0;in;i+) for(j=0;jAij; printf(n输入右端b%d: ,n); for(i=