【最佳平方逼近】hdu4851

来源:互联网 发布:docker nginx 配置 编辑:程序博客网 时间:2024/05/17 09:12


找出不超过n次的多项式f(x)使得上式最小。

设f(x)=sigma(ai*x^i)

将积分式展开有三项,第一项是关于a的二次式,第二项是关于a的一次式,第三项式常数项。

如果想让上式最小,可以对于n+1个系数各求一次偏导,使其偏导为0,则得到n+1个方程,高斯消元即可求解。

关键是这题要用java,但当时还不太会用java,因此解不来。

其实这题是可用希尔波特矩阵直接构造高斯消元,第二天上科学计算课的时候正好就提到了这类题目,推出来和上面求偏导的方法推出来是一样的。

import java.io.*;import java.util.*;import java.math.*;public class Main{    static BigDecimal c[][];    static BigDecimal pro[];    static BigDecimal cal[];    static BigDecimal E;    static int prec=100;    static int r_m;    static int n;//    static BigInteger c[];    private static void swap(int i,int j)    {        BigDecimal tmp;        for (int k=0;k<=n+1;k++) {            tmp=c[i][k];            c[i][k]=c[j][k];            c[j][k]=tmp;        }    }    private static void Gauss()    {        for (int j=0;j<=n;j++) {            int k=j;            BigDecimal Max=new BigDecimal(-1.0);            for (int i=j;i<=n;i++)                if (Max.compareTo(c[i][j].abs())==-1) {                    Max=c[i][j].abs();                    k=i;                }            swap(j,k);            for (int i=0;i<=n;i++)                if (i!=j) {                    BigDecimal tmp;                    tmp=c[i][j].divide(c[j][j],prec,r_m);                    for (k=j+1;k<=n+1;k++)                        c[i][k]=c[i][k].subtract(c[j][k].multiply(tmp));                    c[i][j]=BigDecimal.ZERO;                }        }        for (int i=0;i<=n;i++)            c[i][i]=c[i][n+1].divide(c[i][i],prec,r_m);    }    public static void origin()    {        r_m=BigDecimal.ROUND_DOWN;        pro=new BigDecimal[101];        cal=new BigDecimal[101];        BigDecimal sum=BigDecimal.ZERO;        for (int i=1;i<=100;i++) {            BigDecimal tmp=BigDecimal.valueOf(i);            sum=tmp;            pro[i]=BigDecimal.ONE;            pro[i]=pro[i].divide(sum,prec,r_m);        }        sum=BigDecimal.ONE;        cal[0]=sum;        for (int i=1;i<=100;i++) {            BigDecimal tmp=BigDecimal.valueOf(i);            sum=sum.multiply(tmp);//            cal[i]=BigDecimal.ONE;            cal[i]=sum;        }        E=BigDecimal.ONE;        for (int i=1;i<=100;i++) {            BigDecimal tmp=BigDecimal.ONE;            tmp=tmp.divide(cal[i],prec,r_m);            E=E.add(tmp);        }    }    public static void main(String args[]) throws Exception{        Scanner cin=new Scanner(new BufferedInputStream(System.in));        c=new BigDecimal[100][100];        origin();        while (cin.hasNext()) {            n=cin.nextInt();            for (int i=0;i<=n+1;i++)                for (int j=0;j<=n+1;j++)                    c[i][j]=BigDecimal.ZERO;            for (int i=0;i<=n;i++)                for (int j=0;j<=n;j++) {                        c[i][j]=c[i][j].add(pro[i+j+1]);                        c[j][i]=c[j][i].add(pro[i+j+1]);                }            for (int i=0;i<=n;i++)                 c[i][n+1]=c[i][n+1].subtract(E);            for (int i=0;i<=n;i++)                if (((i+1)&1)!=0) {                    c[i][n+1]=c[i][n+1].add(cal[i]);                }                else {                    c[i][n+1]=c[i][n+1].subtract(cal[i]);                }            for (int i=1;i<=n;i++)                for (int j=i;j<=n;j++)                     if ((i&1)==1) {                        c[j][n+1]=c[j][n+1].add(E.multiply(cal[j].divide(cal[j-i],prec,r_m)));                    }                    else {                        c[j][n+1]=c[j][n+1].subtract(E.multiply(cal[j].divide(cal[j-i],prec,r_m)));                    }            for (int i=0;i<=n;i++)                c[i][n+1]=c[i][n+1].multiply(BigDecimal.valueOf(-2.0));/*        for (int i=0;i<=n;i++) {                for (int j=0;j<=n+1;j++)                    System.out.printf("%.5f ",c[i][j].setScale(55,r_m));                System.out.println();            }*/            Gauss();            for(int i=0;i<=n;i++)               // System.out.println(c[i][i].setScale(60,r_m));                System.out.printf("%.55f\n",c[i][i].setScale(60,r_m));//            System.out.println(pro[3].setScale(55,r_m));                    }    }}


0 0