Given a function f and a set of m>0 distinct points x1 <x2 <⋯<xm . You are supposed to write a function to approximate f by an orthogonal polynomial using the exact function values at the given m points with a weight w(xi ) assigned to each point xi . The total error err=∑i=1m w(xi )[7f(xi )−Pn (xi )]2 must be no larger than a given tolerance.

Format of function:

int OPA( double (*f)(double t), int m, double x[], double w[], double c[], double *eps );

where the function pointer double (*f)(double t) defines the function f; int m is the number of points; double x[] contains points x
1 <x2 <⋯<xm ; double w[] contains the values of a weight function at the given points x[]; double c[] contains the coefficients c0 ,c1 ,⋯,cn of the approximation polynomial Pn (x)=c0 +c
1 x+⋯+cn xn ; double *eps is passed into the function as the tolerance for the error, and is supposed to be returned as the value of error. The function OPA is supposed to return the degree of the approximation polynomial.

Note: a constant Max_n is defined so that if the total error is still not small enough when n = Max_n, simply output the result obtained when n = Max_n.

Sample program of judge:

#include <stdio.h>
#include <math.h>#define MAX_m 200
#define MAX_n 5double f1(double x)
{return sin(x);
}double f2(double x)
{return exp(x);
}int OPA( double (*f)(double t), int m, double x[], double w[], double c[], double *eps );void print_results( int n, double c[], double eps)
{    int i;printf("%d\n", n);for (i=0; i<=n; i++)printf("%12.4e ", c[i]);printf("\n");printf("error = %9.2e\n", eps);printf("\n");
}int main()
{int m, i, n;double x[MAX_m], w[MAX_m], c[MAX_n+1], eps;m = 90;for (i=0; i<m; i++) {x[i] = 3.1415926535897932 * (double)(i+1) / 180.0;w[i] = 1.0;}eps = 0.001;n = OPA(f1, m, x, w, c, &eps);print_results(n, c, eps);m = 200;for (i=0; i<m; i++) {x[i] = 0.01*(double)i;w[i] = 1.0;}eps = 0.001;n = OPA(f2, m, x, w, c, &eps);print_results(n, c, eps);return 0;
}/* Your function will be put here */

Sample Output:
3
-2.5301e-03 1.0287e+00 -7.2279e-02 -1.1287e-01
error = 6.33e-05

4
1.0025e+00 9.6180e-01 6.2900e-01 7.0907e-03 1.1792e-01
error = 1.62e-04

思路

1.需要利用一个结构实现多项式相加，多项式内积。

最简单的方法就是使用数组，并且不需要记录最高次数，只需要将前面的项系数都变为0即可。相当于一个次数为MAX_n的多项式，前面某些项系数为0。

多项式加法：系数直接相加。

多项式乘x:将系数后移一位

多项式乘常数：系数直接乘

多项式内积：比较麻烦的是，变量可能是y也可能是x。一种思路是两个都先按照多项式算出∑k=0max_na[k]⋅xik\sum\limits_{k=0}^{max\_n} a[k]\cdot x_i^kk=0∑max_na[k]⋅xik，之后根据传入的flag确定到底应该选择f(xi)还是∑k=0max_na[k]⋅xik\sum\limits_{k=0}^{max\_n}a[k]\cdot x_i^kk=0∑max_na[k]⋅xik

2.ppt上计算正交多项式的时候用到了ϕ(x)(x−B1)\phi(x)(x-B_1)ϕ(x)(x−B1)，很明显应该拆成ϕ(x)x−B1ϕ(x)\phi(x)x-B_1\phi(x)ϕ(x)x−B1ϕ(x)

3.内积别忘了乘w[i]

4.如果内积没算错，这题大概率就不会错。

上代码

double gx[MAX_m];
double gy[MAX_m];
double gw[MAX_m];
int gm;double Polynomials_Sum(double* f,int flagf, double* g,int flagg)/*计算f(x)，g(y)的求和,flg表示后一个变量*/
{double ret = 0;for (int i = 0; i < gm; i++){   double sumf = 0;double sumg = 0;for (int j = 0; j <= MAX_n; j++){sumf +=  f[j] * pow(gx[i], j);sumg +=  g[j] * pow(gx[i], j);}   sumf = (flagf == 0) ? sumf : gy[i];sumg = (flagg == 0) ? sumg : gy[i];ret += gw[i] * sumg * sumf;}return ret;
}int OPA(double (*f)(double t), int m, double x[], double w[], double c[], double* eps)
{for (int i = 0; i < m; i++){gx[i] = x[i];gw[i] = w[i];gy[i] = f(x[i]);gm = m;}double phi0[MAX_n+1] = { 0 };double phi1[MAX_n+1] = { 0 };double phi2[MAX_n+1] = { 0 };phi0[0] = 1;double y[MAX_n+1] = {0};y[1] = 1;double a0=Polynomials_Sum(phi0, 0,y,1) / Polynomials_Sum(phi0,0, phi0, 0);for (int i = 0; i <= MAX_n; i++){c[i] = phi0[i] * a0;}double error = Polynomials_Sum(y,1, y, 1) - a0 * Polynomials_Sum(phi0,0, y,1 );double xphi0[MAX_n+1] = {0};for (int i = 0; i < MAX_n; i++)xphi0[i+1] = phi0[i];double B1 = Polynomials_Sum(xphi0,0, phi0, 0) / Polynomials_Sum(phi0,0, phi0, 0);phi1[0] = -B1;phi1[1] = 1;a0 = Polynomials_Sum(phi1, 0, y, 1) / Polynomials_Sum(phi1,0, phi1, 0);for (int i = 0; i <= MAX_n; i++){c[i] += a0 * phi1[i];}error -= a0 * Polynomials_Sum(phi1,0, y, 1);int k = 1;while (k < MAX_n && fabs(error) >= *eps){k++;xphi0[0] = 0;for (int i = 0; i < MAX_n; i++)xphi0[i+1] = phi1[i];B1 = Polynomials_Sum(xphi0,0, phi1,0) / Polynomials_Sum(phi1,0, phi1, 0);double C1 = Polynomials_Sum(xphi0,0, phi0, 0) / Polynomials_Sum(phi0,0, phi0, 0);phi2[0] = 0;if (phi1[MAX_n] != 0) {*eps = error;return k;}for (int i = 0; i < MAX_n; i++)phi2[i + 1] = phi1[i];for (int i = 0; i <= MAX_n; i++)phi2[i] = phi2[i] - B1 * phi1[i] - C1 * phi0[i];double t1 = Polynomials_Sum(phi2, 0, y, 1);double t2 = Polynomials_Sum(phi2, 0, phi2, 0);a0 =  t1/t2 ;for (int i = 0; i <= MAX_n; i++)c[i] += a0 * phi2[i];error -= a0 * Polynomials_Sum(phi2,0, y, 1);for (int i = 0; i <= MAX_n; i++){phi0[i] = phi1[i];phi1[i] = phi2[i];}    }*eps = error;return k;
}

6-7 Orthogonal Polynomials Approximation (40分)相关推荐

Java黑皮书课后题第1章：1.12（以千米计的平均速度）假设一个跑步者1小时40分35秒跑了24英里。编写一个程序显示以每小时为多少千米为单位的平均速度值（1英里等于1.6千米）
Java黑皮书课后题第1章:1.12(以千米计的平均速度) 题目题目描述破题代码块修改日志题目题目描述 1.12(以千米计的平均速度)假设一个跑步者1小时40分35秒跑了24英里.编写一个 ...
计算机专业比重点线高40多分,这3所211大学，超过一本线40分就可报考，性价比高，值得报考...
985,211大学是我国的重点大学,与普通大学相比,师资力量和学科力量都很强.所以很多考生之所以在高中期间努力学习,就是为了能考入名校,毕业之后能找一份薪资待遇不错的工作.一般来说,综合实力强的学校录 ...
高中数学40分怎么办_2019年第35届全国高中数学联赛试题及参考答案
2019年第35届全国高中数学联赛考试已结束,本文收集整理本次数学联赛的试题和参考答案,以供大家了解参考. 本次数学联赛由全国高中数学联赛组委会统一命题,共分为一试和二试. 一试时间为80分钟,包括8 ...
高中数学40分怎么办_高二数学不会，准高三该怎么办？40分到高考140如何逆袭？...
原标题:高二数学不会,准高三该怎么办?40分到高考140如何逆袭? 高二,这个年级是有点尴尬的,适应了高一的学习,感觉高二学习没有了动力,离高考还远,于是有些孩子就开始了放任自己,开始了放弃,殊不知高 ...
7-3 素数对猜想 (40 分)
** 7-3 素数对猜想 (40 分) ** 让我们定义d n 为:d n=p n+1 −p n ,其中p i 是第i个素数.显然有d 1 =1,且对于n>1有d n 是偶数."素数 ...
7-2 数组元素循环右移问题 (40 分)
** 7-2 数组元素循环右移问题 (40 分) ** 一个数组A中存有N(>0)个整数,在不允许使用另外数组的前提下,将每个整数循环向右移M(≥0)个位置,即将A中的数据由(A 0 A 1 ⋯ ...
7-1 图形卡片排序游戏 (40 分)
** 7-1 图形卡片排序游戏 (40 分) ** 掌握类的继承.多态性使用方法以及接口的应用.详见[作业指导书2020-OO第07次作业-1指导书V1.0.pdf 输入格式: 首先,在一行上输入一串 ...
7-5 日期问题面向对象设计（聚合二） (40 分)
** 7-5 日期问题面向对象设计(聚合二) (40 分) ** 参考题目7-3的要求,设计如下几个类:DateUtil.Year.Month.Day,其中年.月.日的取值范围依然为:year∈[18 ...
哈登独得40分保罗复出火箭主场103:98复仇魔术
保罗在本场比赛迎来复出. 中新网北京1月28日电北京时间28日,NBA常规赛继续展开争夺,火箭队在主场迎来魔术队的挑战,本场比赛火箭后卫保罗迎来复出.最终凭借哈登独得40分的表现,火箭在主场103: ...

6-7 Orthogonal Polynomials Approximation (40分)

思路

6-7 Orthogonal Polynomials Approximation (40分)相关推荐

最新文章

热门文章