计算机数值分析
实验:
插值方法:
用编程语言编程实现以下算法:
(1) 已知插值节点序列 ,用拉格朗日(Lagrange)插值多项式 计算的函数 在点 的近似值。
(2) 已知插值节点序列 ,用牛顿(Newton)插值值多项式 计算的函数 在点 的近似值。
(3) 用线性函数 拟合给定数据 。
#include<iostream>
#define MAX 100
using namespace std;
//拉格朗日插值
double Lagrange(double x[], double y[], int n, double x0) {
double y0 = 0;
for (int k = 0; k < n; k++) {
double temp = y[k];
for (int i = 0; i < n; i++) {
if (i != k) {
temp =temp * (x0 - x[i]);
temp =temp / (x[k] - x[i]);
}
}
y0 += temp;
}
return y0;
}
//牛顿插值
double Newton(double x[], double y[], int n, double x0) {
//差商表
double temp[MAX][MAX];
for (int i = 1; i < n; i++) {
for (int j = i; j < n; j++) {
if (i == 1) {
temp[i][j] = (y[j] - y[j - 1]) / (x[j] - x[j - 1]);
}
else {
temp[i][j] = (temp[i - 1][j] - temp[i - 1][j - 1]) / (x[j] - x[0]);
}
}
}
double out = 0;
for (int i = 0; i < n; i++) {
if (i == 0)
out += y[0];
else {
double temp1 = temp[i][i];
for (int j = 0; j < i; j++) {
temp1 *= x0 - x[j];
}
out += temp1;
}
}
return out;
}
//最小二乘法拟合数据
void minTwo(double x[], double y[], int n,double &a,double&b) {
double a1=n, b1=0, c1=0;
double a2=0, b2=0, c2=0;
for (int i = 0; i < n; i++) {
b1 += x[i];
c1 += y[i];
a2 += x[i];
b2 += x[i] * x[i];
c2 += x[i] * y[i];
}
a = (c1 * b2 - c2 * b1) / (a1 * b2 - a2 * b1);
b = (c1 * a2 - c2 * a1) / (b1 * a2 - b2 * a1);
}
int main() {
double x[MAX],y[MAX];
int n;
double x0;
cout << "输入插值节点数 n :" << endl;
cin >> n;
for (int i = 0; i < n; i++) {
cout << i+1 << ":";
cin >> x[i] >> y[i];
}
//cout << "输入要计算的 x0 :";
//cin >> x0;
//cout << "函数f(x)的近似值为:Lagrange:" << Lagrange(x, y, n, x0) << endl;
//cout << "函数f(x)的近似值为:Newton:" << Newton(x, y, n, x0) << endl;
double a, b;
minTwo(x, y, n,a,b);
cout << "拟合函数为:y = " << a << " + " << b << "x " << endl;
return 0;
}
数值积分与微分
用编程语言编程实现以下算法:
(1) 用复化梯形公式 的自动控制误差算法求积分 。
(2) Romberg积分算法求积分 。
#include<iostream>
#define MAX 20
using namespace std;
//测试函数
double f(double x) {
/*if (x == 0)
return 1;
double a = sin(x);
a = a / x;*/
return sqrt(x*x*x);
}
//复合梯形公式
double Trapz(double a,double b,int n,double(*f)(double x)) {
double h = (b - a) / n;
double out = 0;
out += f(a)/2;
out += f(b)/2;
for (int i = 1; i < n; i++) {
out += f(a + i * h);
}
out *= h;
return out;
}
//Remberg算法
double Remberg(double a,double b,int n,double (*f)(double x)) {
//记录数据
double temp[MAX][MAX] = {0};
temp[0][0] = Trapz(a, b, 1, f);
for (int k = 1; k <= n; k++) {
for (int m = 0; m <= k; m++) {
if (m == 0) {
double hLast = (b - a) / pow(2, k - 1);
double p = 0;
for (int i = 0; i < pow(2, k - 1); i++) {
p += f(a + i * hLast + 0.5 * hLast);
}
p = p * hLast / 2;
temp[k][m] = temp[k - 1][m] / 2 + p;
}
else {
temp[k][m] = (pow(4, m) * temp[k][m - 1] - temp[k - 1][m - 1]) / (pow(4, m) - 1);
}
cout << temp[k][m]<<"\t";
}
cout << endl;
}
return temp[n][n];
}
int main() {
//测试数据
//cout << Trapz(0, 1, 1, f) << endl;
cout << Remberg(0, 1, 5, f) << endl;
return 0;
}
本章学习了数值积分的算法,在计算一些无法通过直接计算来求解的函数积分时,来实现各种近似积分的方法,包含插值型求积公式、Newton-Cotes公式、复合梯形公式、复合Simpson公式,以及Romberg算法等。在运用计算机求解积分时的技巧和方法。在掌握这些东西的时候虽然难以理解,但这些东西都包含数学的魅力。
常微分方程初值问题的数值解法
用编程语言编程实现以下算法:
(1) 用改进的欧拉(Euler)公式求解常微分方程初值问题。
(2) 用四阶龙格-库塔(Runge-Kutta)方法求解常微分方程初值问题。
#include<iostream>
using namespace std;
//导数
double f_xy(double x, double y) {
//return -1 * y + x + 1;
return -1*y;
}
//欧拉算法
double Euler(double x0,double y0, double h, double xn,double(*f)(double x,double y)) {
int n = (xn - x0) / h;
double yNext=y0;
double xNext = x0;
for (int i = 0; i < n; i++) {
yNext = yNext + h / 2 * (f_xy(xNext, yNext) + f(xNext + h, yNext + h * f(xNext, yNext)));
xNext += h;
}
return yNext;
}
//龙格库塔算法
double RungeKutta(double x0,double y0,double h,double xn,double(*test)(double x,double y)) {
int n = (xn - x0) / h;
double xNext = x0;
double yNext = y0;
for (int i = 0; i < n; i++) {
double k1 = f_xy(xNext, yNext);
double k2 = f_xy(xNext + h / 2, yNext + h / 2 * k1);
double k3 = f_xy(xNext + h / 2, yNext + h / 2 * k2);
double k4 = f_xy(xNext + h, yNext + h * k3);
xNext += h;
yNext += h / 6 * (k1 + 2 * k2 + 2 * k3 + k4);
}
return yNext;
}
int main() {
cout << "改进的欧拉算法:" << endl;
cout << Euler(0, 1,0.05, 0.1, f_xy) << endl;
cout << Euler(0, 1,0.05, 0.2, f_xy) << endl;
cout << Euler(0, 1,0.05, 0.3, f_xy) << endl;
cout << Euler(0, 1,0.05, 0.4, f_xy) << endl;
cout << Euler(0, 1,0.05, 0.5, f_xy) << endl;
cout << "龙格库塔算法:" << endl;
cout << RungeKutta(0, 1, 0.1, 0.1, f_xy) << endl;
cout << RungeKutta(0, 1, 0.1, 0.2, f_xy) << endl;
cout << RungeKutta(0, 1, 0.1, 0.3, f_xy) << endl;
cout << RungeKutta(0, 1, 0.1, 0.4, f_xy) << endl;
cout << RungeKutta(0, 1, 0.1, 0.5, f_xy) << endl;
return 0;
}
本章学习的是微分方程的近似解法,已知初值(x0,y0),和求积函数,求解xn时yn的近似值,问题是对于不能直接通过计算得出结果计算式子。总体思想就是求曲线的积分,求面积来代替求解,然后使用微分法,离散化成一个一个区域求面积,所有面积的和就是yn,每个面积用长方形代替长取区域长度,宽取区域左端点求的方法交欧拉算法(Euler),取右端点的叫隐式欧拉算法,取中点的叫中点欧拉格式,用梯形计算区域面积的叫梯形格式,由于右端点未知,用欧拉算法先求出来,再用梯形公式校正,就是改进的欧拉算法。龙格库塔算法就是在区域内取几个点,把对应的y加权求平均作为长方形的高。求解。
方程求根的数值方法
用编程语言编程实现以下算法:
(1) 用二分法求 的根。
(2) 用牛顿(Newton)迭代法求 在 附近的根。
(3) 用弦截法求 的根。
#include<iostream>
#define MAX 1000
using namespace std;
//求解函数
double f(double x) {
return x*x*x-x-1;
}
//导函数
double fl(double x) {
return 2*x+cos(x);
}
//二分法
double dichotomy(double a,double b,double precision,double(*f)(double x)) {
if (f(a) * f(b) > 0)
return NULL;
else {
while (1) {
double temp = (a + b) / 2;
if (b - a < precision){
return temp;
break;
}
else if (f(a) * f(temp) < 0) {
b = temp;
}
else {
a = temp;
}
}
}
}
//Newton迭代法
double Newton(double x0, double precision,double(*f)(double x),double(*fl)(double x)) {
double xn = x0;
for (int i = 0; i < MAX; i++) {
xn = xn - f(xn) / fl(xn);
if (f(xn) < precision&&-f(xn)<precision)
return xn;
}
return NULL;
}
//弦切法
double stringCut(double x0,double x1, double precision, double(*f)(double x)) {
double xn;
for (int i = 0; i < MAX; i++) {
xn = x1 - f(x1) / (f(x1)-f(x0))*(x1-x0);
if (xn-x1 < precision && x1-xn < precision)
return xn;
x0 = x1;
x1 = xn;
}
return NULL;
}
int main() {
//cout << dichotomy(1, 2, 0.01, f) << endl;
//cout<<Newton(0.5, 0.001, f, fl) << endl;
//cout<<Newton(-1.5, 0.001, f, fl) << endl;
cout <<"结果:"<< stringCut(1,2,0.00001,f) << endl;
return 0;
}