计算机数值分析-实验-用C++实现

计算机数值分析

实验:

实验报告
代码文件

插值方法:

用编程语言编程实现以下算法:
(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;
}
  • 6
    点赞
  • 37
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值