数值计算基础

ゞ浴缸里的玫瑰 2022-06-07 10:49 397阅读 0赞

**1:第二章**

**(1)二分法求解非线性方程：**

\#include<stdio.h>

\#include <math.h>

\#define f(x) ((x\*x-1)\*x-1)

void main()

\{ float a,b,x,eps;

int k=0;

printf("intput eps\\n");/\* 容许误差 \*/

scanf("%f",&eps);

printf("a,b=\\n");

for(;;)

\{scanf("%f, %f",&a ,&b);

if(f(a)\*f(b)>=0) /\* 判断是否符合二分法使用的条件 \*/

printf("二分法不可使用,请重新输入:\\n");

else break;

\{ x=(a+b)/2;

k++;

if(f(a)\*f(x)<0) /\* 如果f(a)\*f(x)<0，则根在区间的左半部分 \*/

b=x;

else if(f(a)\*f(x)>0) /\* 否则根在区间的右半部分 \*/

a=x;

else break;

\}while(fabs(b-a)>eps);/\*判断是否达到精度要求,若没有达到,继续循环\*/

x=(a+b)/2; /\* 取最后的小区间中点作为根的近似值 \*/

printf("\\n The root is x=%f, k=%d\\n",x,k);

**运行结果：**

intput eps

0.00001

a,b=

2,-5

Theroot is x=1.324721, k=20

Press any key to continue

总结：本题关键在于两个端点的取值和误差的判断，此程序较容易。二分法收敛速度较快，但缺点是只能求解单根。

**(2)牛顿法求解非线性方程：**

\#include <stdio.h>

\#include <math.h>

float f(float x) /\* 定义函数f(x) \*/

\{ return((-3\*x+4)\*x-5)\*x+6; \}

float f1(float x) /\* 定义函数f(x)的导数 \*/

\{ return (-9\*x+8)\*x-5; \}

void main()

\{ float eps,x0,x1=1.0;

printf("input eps:\\n");

scanf("%f",&eps); /\* 输入容许误差 \*/

\{ x0=x1; /\* 准备下一次迭代的初值 \*/

x1=x0-f(x0)/f1(x0); /\* 牛顿迭代 \*/

\}while(fabs(x1-x0)>eps); /\*当满足精度,输出近似根\*/

printf("x=%f\\n",x1);

**程序运行结果：**

x=1.265328

总结：关键是牛顿迭代的应用，程序中最大缺点是函数及其导数已唯一给出确定不可求的随意函数的根，牛顿法比二分法收敛快，可以求重根。

**2:第三章**

**(1)列主元素消去法求解线性方程：**

\#include<iostream>

\#include<cmath>

\#define N 20

using namespace std;

void load();

float a\[N\]\[N\];

int m;

int main()\{

inti,j;

intc,k,n,p,r;

floatx\[N\],l\[N\]\[N\],s,d;

cout<<"下面请输入未知数的个数m=";

cin>>m;

cout<<endl;

cout<<"请按顺序输入增广矩阵a:"<<endl;

load();

for(i=0;i<m;i++)

for(j=i;j<m;j++)

c=(fabs(a\[j\]\[i\])>fabs(a\[i\]\[i\]))?j:i; /\*找列最大元素\*/

for(n=0;n<m+1;n++)

\{s=a\[i\]\[n\]; a\[i\]\[n\]=a\[c\]\[n\]; a\[c\]\[n\]=s;\} /\*将列最大数防在对角线上\*/

for(p=0;p<m+1;p++)

cout<<a\[i\]\[p\]<<"\\t";

cout<<endl;

for(k=i+1;k<m;k++)

l\[k\]\[i\]=a\[k\]\[i\]/a\[i\]\[i\];

for(r=i;r<m+1;r++) /\*化成三角阵\*/

a\[k\]\[r\]=a\[k\]\[r\]-l\[k\]\[i\]\*a\[i\]\[r\];

x\[m-1\]=a\[m-1\]\[m\]/a\[m-1\]\[m-1\];

for(i=m-2;i>=0;i--)

d=0;

for(j=i+1;j<m;j++)

d=d+a\[i\]\[j\]\*x\[j\];

x\[i\]=(a\[i\]\[m\]-d)/a\[i\]\[i\]; /\*求解\*/

cout<<"该方程组的解为:"<<endl;

for(i=0;i<m;i++)

cout<<"x\["<<i<<"\]="<<x\[i\]<<"\\t";

//system("pause");

return0;

void load()

int i,j;

for(i=0;i<m;i++)

for(j=0;j<m+1;j++)

cin>>a\[i\]\[j\];

**运行结果：**

下面请输入未知数的个数m=3

请按顺序输入增广矩阵a:

1 2 3 4

5 1 0 8

4 6 9 2

0 -6.5 -11.25 5.5

0 -1.86265e-008 -0.115385 3.92308

该方程组的解为:

x\[0\]=-9.99999 x\[1\]=58 x\[2\]=-34 Press any key to continue

总结：列主元素消去法的目的是为了防止减去一个较小的数时大数淹没小数，而使结果产生较大误差，本程序关键在每次消元时找到相应列中的最大项，然后交换两行位置，在进行计算。

**(2)LU分解法求解线性方程：**

\#include<stdio.h>

voidsolve(float l\[\]\[100\],float u\[\]\[100\],float b\[\],float x\[\],int n)

\{inti,j;

floatt,s1,s2;

floaty\[100\];

for(i=1;i<=n;i++)/\* 第一次回代过程开始 \*/

\{s1=0;

for(j=1;j<i;j++)

t=-l\[i\]\[j\];

s1=s1+t\*y\[j\];

y\[i\]=(b\[i\]+s1)/l\[i\]\[i\]; \}

for(i=n;i>=1;i--)/\* 第二次回代过程开始 \*/

s2=0;

for(j=n;j>i;j--)

t=-u\[i\]\[j\];

s2=s2+t\*x\[j\];

x\[i\]=(y\[i\]+s2)/u\[i\]\[i\];

voidmain()

\{floata\[100\]\[100\],l\[100\]\[100\],u\[100\]\[100\],x\[100\],b\[100\];

inti,j,n,r,k;

floats1,s2;

for(i=1;i<=99;i++)/\*将所有的数组置零，同时将L矩阵的对角值设为1\*/

for(j=1;j<=99;j++)

l\[i\]\[j\]=0,u\[i\]\[j\]=0;

if(j==i) l\[i\]\[j\]=1;

printf("input n:\\n");/\*输入方程组的个数\*/

scanf("%d",&n);

printf("input array A:\\n");/\*读取原矩阵A\*/

for(i=1;i<=n;i++)

for(j=1;j<=n;j++)

scanf("%f",&a\[i\]\[j\]);

printf("input array B:\\n");/\*读取列矩阵B\*/

for(i=1;i<=n;i++)

scanf("%f",&b\[i\]);

for(r=1;r<=n;r++)/\*求解矩阵L和U\*/

for(i=r;i<=n;i++)

s1=0;

for(k=1;k<=r-1;k++)

s1=s1+l\[r\]\[k\]\*u\[k\]\[i\];

u\[r\]\[i\]=a\[r\]\[i\]-s1;

for(i=r+1;i<=n;i++)

\{s2=0;

for(k=1;k<=r-1;k++)

s2=s2+l\[i\]\[k\]\*u\[k\]\[r\];

l\[i\]\[r\]=(a\[i\]\[r\]-s2)/u\[r\]\[r\];

printf("arrayL:\\n");/\*输出矩阵L\*/

for(i=1;i<=n;i++)

for(j=1;j<=n;j++)

printf("%7.3f ",l\[i\]\[j\]);

printf("\\n");

printf("arrayU:\\n");/\*输出矩阵U\*/

for(i=1;i<=n;i++)

for(j=1;j<=n;j++)

printf("%7.3f ",u\[i\]\[j\]);

printf("\\n");

solve(l,u,b,x,n);

printf("解为:\\n");

for(i=1;i<=n;i++)

printf("x%d=%f\\n",i,x\[i\]);

**运行结果：**

input n:

inputarray A:

2 2 3

4 7 7

\-2 4 5

inputarray B:

3 1 -7

array L:

1.000 0.000 0.000

2.000 1.000 0.000

\-1.000 2.000 1.000

array U:

2.000 2.000 3.000

0.000 3.000 1.000

0.000 0.000 6.000

解为:

x1=2.000000

x2=-2.000000

x3=1.000000

Pressany key to continue

总结：关键是把矩阵分解为L、U两个三角矩阵，回代过程比较简单。

**3:第四章**

**(1)拉格朗日差值多项式;**

\#include<stdio.h>

\#include<math.h>

\#define MAX 100

void main()

\{ int i,j,k,m,n,N,mi;

float tmp,mx;

float X\[MAX\]\[MAX\],Y\[MAX\],x\[MAX\],y\[MAX\],a\[MAX\];

printf("\\n 输入拟合多项式的次数:\\n");

scanf("%d",&m);

printf("\\n 输入给定点的个数n及坐标(x,y):\\n");

scanf("%d",&N);

printf("\\n");

for(i=0;i<N;i++)

scanf("%f,%f",&x\[i\],&y\[i\]);

for(i=0;i<=m;i++)

for(j=i;j<=m;j++)

tmp=0;

for(k=0;k<N;k++)

tmp=tmp+pow(x\[k\],(i+j));

X\[i\]\[j\]=tmp;

X\[j\]\[i\]=X\[i\]\[j\];

for(i=0;i<=m;i++)

tmp=0;

for(k=0;k<N;k++)

tmp=tmp+y\[k\]\*pow(x\[k\],i);

Y\[i\]=tmp;

for(j=0;j<m;j++)

for(i=j+1,mi=j,mx=fabs(X\[j\]\[j\]);i<=m;i++)

if(fabs(X\[i\]\[j\])>mx)

mi=i;

mx=fabs(X\[i\]\[j\]);

if(j<mi)

tmp=Y\[j\];

Y\[j\]=Y\[mi\];

Y\[mi\]=tmp;

for(k=j;k<=m;k++)

tmp=X\[j\]\[k\];

X\[j\]\[k\]=X\[mi\]\[k\];

X\[mi\]\[k\]=tmp;

for(i=j+1;i<=m;i++)

tmp=-X\[i\]\[j\]/X\[j\]\[j\];

Y\[i\]+=Y\[j\]\*tmp;

for(k=j;k<=m;k++)

X\[i\]\[k\]+=X\[j\]\[k\]\*tmp;

a\[m\]=Y\[m\]/X\[m\]\[m\];

for(i=m-1;i>=0;i--)

a\[i\]=Y\[i\];

for(j=i+1;j<=m;j++)

a\[i\]-=X\[i\]\[j\]\*a\[j\];

a\[i\]/=X\[i\]\[i\];

printf("\\n所求的二次多项式为:\\n");

printf("P(x)=%f",a\[0\]);

for(i=1;i<=m;i++)

printf("+(%f)\*x^%d",a\[i\],i);

**运行结果：**

输入拟合多项式的次数:

输入给定点的个数n及坐标(x,y):

1,2

5,3

4,2

所求的二次多项式为:

P(x)=1.980417+(0.282759)\*x^1+(-0.299937)\*x^2+(0.022071)\*x^3+(0.016624)\*x^4+(-0.0

01934)\*x^5Press any key to continue

总结：拉格朗日计算公式中，只需要知道各个点即可

**4:第五章**

**(1)曲线拟合：**

\#include<stdio.h>

\#include<math.h>

\#define MAX 100

void main()

\{ int i,j,k,m,n,N,mi;

float tmp,mx;

float X\[MAX\]\[MAX\],Y\[MAX\],x\[MAX\],y\[MAX\],a\[MAX\];

printf("\\n 输入拟合多项式的次数:\\n");

scanf("%d",&m);

printf("\\n 输入给定点的个数n及坐标(x,y):\\n");

scanf("%d",&N);

printf("\\n");

for(i=0;i<N;i++)

scanf("%f,%f",&x\[i\],&y\[i\]);

for(i=0;i<=m;i++)

for(j=i;j<=m;j++)

tmp=0;

for(k=0;k<N;k++)

tmp=tmp+pow(x\[k\],(i+j));

X\[i\]\[j\]=tmp;

X\[j\]\[i\]=X\[i\]\[j\];

for(i=0;i<=m;i++)

tmp=0;

for(k=0;k<N;k++)

tmp=tmp+y\[k\]\*pow(x\[k\],i);

Y\[i\]=tmp;

for(j=0;j<m;j++)

for(i=j+1,mi=j,mx=fabs(X\[j\]\[j\]);i<=m;i++)

if(fabs(X\[i\]\[j\])>mx)

mi=i;

mx=fabs(X\[i\]\[j\]);

if(j<mi)

tmp=Y\[j\];

Y\[j\]=Y\[mi\];

Y\[mi\]=tmp;

for(k=j;k<=m;k++)

tmp=X\[j\]\[k\];

X\[j\]\[k\]=X\[mi\]\[k\];

X\[mi\]\[k\]=tmp;

for(i=j+1;i<=m;i++)

tmp=-X\[i\]\[j\]/X\[j\]\[j\];

Y\[i\]+=Y\[j\]\*tmp;

for(k=j;k<=m;k++)

X\[i\]\[k\]+=X\[j\]\[k\]\*tmp;

a\[m\]=Y\[m\]/X\[m\]\[m\];

for(i=m-1;i>=0;i--)

a\[i\]=Y\[i\];

for(j=i+1;j<=m;j++)

a\[i\]-=X\[i\]\[j\]\*a\[j\];

a\[i\]/=X\[i\]\[i\];

printf("\\n所求的二次多项式为:\\n");

printf("P(x)=%f",a\[0\]);

for(i=1;i<=m;i++)

printf("+(%f)\*x^%d",a\[i\],i);

**输入拟合多项式的次数:**

**2**

**输入给定点的个数n及坐标(x,y):**

**5**

**1,2**

**5,3**

**2,4**

**8,3**

**-1,5**

**所求的二次多项式为:**

**P(x)=3.952280+(-0.506315)\*x^1+(0.050877)\*x^2Pressany key to continue**

**5:第六章**

**(1)辛普生求积方法：**

\#include <stdio.h>

\#define N 16 /\* 等分数 \*/

float func(float x)

\{ float y;

y=4.0/(1+x\*x);

return(y);

void gedianzhi(float y\[\],float a,float h)

\{ int i;

for(i=0;i<=N;i++)

y\[i\]=func(a+i\*h);

float simpson(float y\[\],float h)

\{ float s,s1,s2;

int i;

s1=y\[1\];

s2=0.0;

for(i=2;i<=N-2;i=i+2)

\{ s1+=y\[i+1\]; /\* 计算奇数项的函数值之和 \*/

s2+=y\[i\]; /\* 计算偶数项的函数值之和 \*/

s=y\[0\]+y\[N\]+4.0\*s1+2.0\*s2;

return(s\*h/3.0);

main()

\{ float a,b,h,s,f\[N+1\];

scanf("%f,%f",&a,&b);

h=(b-a)/( float)N;

gedianzhi(f,a,h);

s=simpson(f,h);

printf("s=%f\\n",s);

**运行结果：**

1,3

s=1.854590

Press any key to continue

总结：辛普生算法是一种积分方法，采用三点法插值，如果h较小的话，误差很小，因为它的插值余项，辛普生算法比较精确，程序关键是对所取的点的取和，注意

**6:第七章**

(1)改进欧拉法求解常微分方程的初值问题

\#include <stdio.h>

float func(float x,float y)

\{ return(y-x);

float euler(float x0,float xn,float y0,intN)

\{ float x,y,yp,yc,h;

int i;

x=x0;

y=y0;

h=(xn-x0)/(float)N;

for(i=1;i<=N;i++)

\{ yp=y+h\*func(x,y);

x=x0+i\*h;

yc=y+h\*func(x,yp);

y=(yp+yc)/2.0;

return(y);

main()

\{ float x0,xn,y0,e;

int n;

printf("\\ninput n:\\n ");

scanf("%d",&n);

printf("input x0,xn:\\n ");

scanf("%f,%f",&x0,&xn);

printf("input y0:\\n ");

scanf("%f",&y0);

e=euler(x0,xn,y0,n);

printf("y(%f)=%6.4f",y0,e);

**inputn:**

**20**

**inputx0,xn:**

**1,6**

**inputy0:**

**2**

**y(2.000000)=7.0000Pressany key to continue**

(2)四阶龙格—库塔法

\#include <stdio.h>

float func(float x,float y)

\{ return(x-y);

float runge\_kutta(float x0,float xn,floaty0,int N)

\{ float x,y,y1,y2,h,xh;

float d1,d2,d3,d4;

int i;

x=x0;

y=y0;

h=(xn-x0)/(float)N;

for(i=1;i<=N;i++)

\{ xh=x+h/2;

d1=func(x,y);

d2=func(xh,y+h\*d1/2.0);

d3=func(xh,y+h\*d2/2.0);

d4=func(xh,y+h\*d3);

y=y+h\*(d1+2\*d2+2\*d3+d4)/6.0;

x=x0+i\*h; \}

return(y);

main()

\{ float x0,xn,y0,e;

int N;

printf("\\ninput n:\\n ");

scanf("%d",&N);

printf("input x0,xn:\\n ");

scanf("%f,%f",&x0,&xn);

printf("input y0:\\n ");

scanf("%f",&y0);

e=runge\_kutta(x0,xn,y0,N);

printf("y(%f)=%8.6f",y0,e);

**input n:**

**10**

**input x0,xn:**

**1,2**

**input y0:**

**5**

**y(5.000000)=2.833863Press any key to continue**

**2-2Gauss-Seidel方法**

\#include<math.h>

\#include<stdio.h>

intgsdl(a,b,n,x,eps)

int n;

doublea\[\],b\[\],x\[\],eps;

\{int i,j,u,v;

double p,t,s,q;

for(i=0;i<=n-1;i++)

\{u=i\*n+i;

p=0.0;

x\[i\]=0.0;

for(j=0;j<=n-1;j++)

\{v=i\*n+j;

p=p+fabs(a\[v\]);

if(p>=fabs(a\[u\]))

\{printf(“fail\\n”);

return(-1);

p=eps+1.0;

while(p>=eps)

\{for(i=0;i<=n-1;i++)

\{t=x\[i\];s=0.0;

for(j=0;j<=n-1;j++)

\{if(j!=i)

\{s=s+a\[i\*n+j\]\*x\[j\];

x\[i\]=(b\[i\]-s)/a\[i\*n+j\];

q=fabs(x\[i\]-t)/(1.0+fabs(x\[i\]));

if(q>p)

\{p=q;

return(1);

main()

\{int i;

double eps;

static doublea\[4\]\[4\]=\{ \{7,2,1,-2\}\{9,15,3,-2\}\{-2,-2,11,5\}\{1,3,2,13\}\};

static doublex\[5\],b\[4\]=\{4,7,-1,0\};

eps=0.000001;

if(dsdl(a,b,4,x,eps)>0)

\{for(i=0;i<=3;i++)

\{printf(“x(%d)=%13.7e\\n”,i,x\[i\]);

数值计算基础

发表评论取消回复

还没有评论，来说两句吧...

相关阅读

相关 python学习——NumPy数值计算基础

相关 Shell 数值计算

相关 MATLAB07:数值计算

相关 java数值计算_Java数值计算包 Commons Math

相关利用计算机进行复杂问题的数值计算,深度学习之数学基础（数值计算）

相关【基础】数值类型

相关数值计算——多项式插值

相关数值计算——求解非线性方程组

相关数值计算基础

相关 Scipy——数值计算库