數值計算基礎
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;
}
do
{ 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); /* 輸入容許誤差 */
do
{ 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
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:
3
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);
}
執行結果：
輸入擬合多項式的次數:
5
輸入給定點的個數n及座標(x,y):
3
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]);
}
}
相關文章
- MapReduce程式設計基礎（二）——數值概要（計算最大值、最小值、平均值）程式設計
- Octave 數值計算
- 雲端計算基礎
- 計算機基礎計算機
- 數值計算 插值與擬合
- Python基礎學習篇-2-數值運算和字串Python字串
- Python3 基礎學習之數值簡單運算Python
- python計算對數值Python
- 計算機基礎：位運算計算機
- 計算機基礎-網路基礎計算機
- 計算機基礎-Socket計算機
- 大學計算機基礎計算機
- 雲端計算基礎-0
- Java 基礎 之 算數運算子Java
- javascript基礎（算數運算子）（八）JavaScript
- 基礎1：Python 數值型別Python型別
- 計算機基礎:離散數學和完備性計算機
- 計算機系統002 – 數值運算計算機
- 計算機系統002 - 數值運算計算機
- 【scipy 基礎】--空間計算
- 計算機基礎系列一：計算機硬體計算機
- 數值計算的可靠性（一）
- 數值計算的可靠性（二）
- 數值計算的可靠性（三）
- 【ES6基礎】預設引數值
- Go基礎系列：惰性數值生成器Go
- 05-JavaScript基礎-算數運算子JavaScript
- TensorFlow 計算與智慧基礎
- 計算機網路---基礎篇計算機網路
- Shell程式設計-04-Shell中變數數值計算程式設計變數
- 計算機基礎之計算機硬體系統計算機
- 1.2程式設計基礎之變數定義、賦值及轉換程式設計變數賦值
- 【C基礎】整形提升與算數轉換
- 數理統計基礎 統計量
- Linux基礎命令---ipcalc計算iPLinuxPCA
- 學計算機需要什麼基礎計算機
- 學計算機需要什麼基礎？計算機
- 基礎面試題 — 計算機網路面試題計算機網路