数值分析程序

合集下载

数值分析插值算法源程序

数值分析插值算法源程序

#include<stdio.h>#include<math.h>float f(float x) //计算ex的值{return (exp(x));}float g(float x) //计算根号x的值{return (pow(x,0.5));}void linerity () //线性插值{float px,x;float x0,x1;printf("请输入x0,x1的值\n");scanf("%f,%f",&x0,&x1);printf("请输入x的值: ");scanf("%f",&x);px=(x-x1)/(x0-x1)*f(x0)+(x-x0)/(x1-x0)*f(x1);printf("f(%f)=%f \n",x,px);}void second () //二次插值{float x0,x1,x2,x,px;x0=0;x1=0.5;x2=2;printf("请输入x的值:");scanf("%f",&x);px=((x-x1)*(x-x2))/((x0-x1)*(x0-x2))*f(x0)+((x-x0)*(x-x2))/((x1-x0)*(x1-x2))*f(x1)+((x-x0)* (x-x1))/((x2-x0)*(x2-x1))*f(x2);printf("f(%f)=%f\n",x,px);}void Hermite () //Hermite插值{int i,k,n=2;int flag1=0;printf("Hermite插值多项式H5(x)=");for(i=0;i<=n;i++){int flag=0;flag1++;if(flag1==1){printf("y%d[1-2(x-x%d)*(",i,i);}else{printf("+y%d[1-2(x-x%d)*(",i,i);}for(k=0;k<=n;k++){if(k!=i){flag++;if(flag==1){printf("(1/x%d-x%d)",i,k);}else{printf("+(1/x%d-x%d)",i,k);}}}printf(")]*(");for(k=0;k<=n;k++){if(i!=k){printf("[(x-x%d)/(x%d-x%d)]2",i,k,i);}}printf(")");}printf("\n");}void sectionl () //分段线性插值{float x[5]={2.0,2.1,2.2,2.3,2.4};float y;printf("请输入y:");scanf("%f",&y);if(y>=2.0&&y<2.1){float px;px=((y-x[1])/(x[0]-x[1]))*g (x[0])+((y-x[0])/(x[1]-x[0]))*g (x[1]);printf("f(%f)=%f\n",y,px);}else if(y>=2.1&&y<2.2){float px;px=((y-x[2])/(x[1]-x[2]))*g (x[1])+((y-x[1])/(x[2]-x[1]))*g (x[2]);printf("f(%f)=%f\n",y,px);}else if(y>=2.2&&y<2.3){float px;px=((y-x[3])/(x[2]-x[3]))*g (x[2])+((y-x[2])/(x[3]-x[2]))*g (x[3]);printf("f(%f)=%f\n",y,px);}else if(y>=2.3&&y<2.4){float px;px=((y-x[4])/(x[3]-x[4]))*g (x[3])+((y-x[3])/(x[4]-x[3]))*g (x[4]);printf("f(%f)=%f\n",y,px);}else if(y>2.4) printf("**********ERROR!******************\n"); }void sectionp (){int i;float a[5]={2.0,2.1,2.2,2.3,2.4};float x,y;printf("input the data: x?\n");scanf("%f",&x);if(x<a[1]){i=1;goto loop;}if(x>a[4]){i=4;goto loop;}i=1;loop1:i++;if(x>a[i])goto loop1;if(fabs(x-a[i-1])<=fabs(x-a[i]))i=i-1;loop:y=g(a[i-1])*(x-a[i])*(x-a[i+1])/((a[i-1]-a[i])*(a[i-1]-a[i+1]));y=y+g(a[i])*(x-a[i-1])*(x-a[i+1])/((a[i]-a[i-1])*(a[i]-a[i+1]));y=y+g(a[i+1])*(x-a[i-1])*(x-a[i])/((a[i+1]-a[i-1])*(a[i+1]-a[i]));printf("f(%f)=%f\n",x,y);}int main(){char flag1='y';while(flag1=='y'){int flag=0;printf("*******[1]:线性插值***************\n");printf("*******[2]:二次插值***************\n");printf("*******[3]:Hermite插值************\n");printf("*******[4]:分段线性插值***********\n");printf("*******[5]:分段抛物线插值*********\n");printf("请输入:");scanf("%d",&flag);switch(flag){case 1:linerity ();break;case 2:second ();break;case 3:Hermite ();break;case 4:sectionl ();break;case 5:sectionp ();break;default:printf("error!!\n");}printf("是否继续?y/n \n");getchar();scanf("%c",&flag1);}return 0;}。

数值分析方法

数值分析方法

数值分析方法数值分析方法是一种通过数学模型和计算方法来解决实际问题的技术。

它在科学计算、工程设计、经济分析等领域有着广泛的应用。

数值分析方法的核心在于将连续的数学问题转化为离散的计算问题,通过数值计算来逼近解析解,从而得到问题的近似解。

本文将介绍数值分析方法的基本原理、常用技术和应用领域。

数值分析方法的基本原理是利用数值计算来逼近解析解。

在实际问题中,很多数学模型很难或者无法得到精确的解析解,这时就需要借助数值分析方法来求解。

数值分析方法的基本步骤包括建立数学模型、离散化、选择适当的数值计算方法、计算近似解并进行误差分析。

其中,离散化是数值分析方法的核心,它将连续的数学问题转化为离散的计算问题,从而使得问题可以通过计算机进行求解。

常用的数值分析方法包括插值法、数值积分、常微分方程数值解、偏微分方程数值解等。

插值法是一种通过已知数据点来估计未知数据点的方法,常用的插值方法包括拉格朗日插值、牛顿插值等。

数值积分是一种通过数值计算来逼近定积分的方法,常用的数值积分方法包括梯形法则、辛普森法则等。

常微分方程数值解和偏微分方程数值解是解决微分方程数值解的常用方法,常用的数值解方法包括欧拉法、龙格-库塔法等。

数值分析方法在科学计算、工程设计、经济分析等领域有着广泛的应用。

在科学计算中,数值分析方法常用于模拟物理现象、计算数学模型等。

在工程设计中,数值分析方法常用于求解结构力学、流体力学等问题。

在经济分析中,数值分析方法常用于求解经济模型、金融衍生品定价等问题。

总之,数值分析方法已经成为现代科学技术和工程技术中不可或缺的一部分。

综上所述,数值分析方法是一种通过数学模型和计算方法来解决实际问题的技术。

它的基本原理是利用数值计算来逼近解析解,常用的方法包括插值法、数值积分、常微分方程数值解、偏微分方程数值解等。

数值分析方法在科学计算、工程设计、经济分析等领域有着广泛的应用。

希望本文的介绍能够帮助读者更好地理解数值分析方法的基本原理和应用价值。

数值分析计算方法程序汇总

数值分析计算方法程序汇总

(一)复化梯形公式例:求121?x dx -=⎰程序:#include "stdio.h"void main(){double a,b,s,h,x;int i,n;a=-1.0;b=1.0;n=10;h=(b-a)/n;x=a;s=x*x/2;for(i=1;i<n;i++){x=x+h;s=s+x*x;}s=s+b*b/2;s=s*h;printf("s=%f\n",s);}结果:s=0.680000(二)复化辛普森公式例:求130?x dx=⎰程序:#include "stdio.h"void main(){double a,b,c,s,h,x,y;int i,n;a=0.0;b=1.0;n=10;s=0.0;h=(b-a)/n;x=a;y=x+h;c=(x+y)/2;for(i=1;i<=n;i++){s=s+x*x*x+4*c*c*c+y*y*y;x=x+h;y=y+h;c=c+h;}s=s*h/6;printf("s=%f\n",s);}结果:s=0.250000(三)复化高斯公式例:求220?x dx=⎰程序:#include <stdio.h>#include <math.h>main(){double a,b,h,s,x1,x2;int i,n;a=0;b=2;n=20;s=0;h=(b-a)/n;for(i=0;i<n;i++){x1=a+i*h+h/2*(1/1.732+1); x2=a+i*h+h/2*(1-1/1.732); s=s+x1*x1*x1+x2*x2*x2; }s=h/2*s;printf("s=%f\n",s);}结果:s=4.000000(四)迭代法例:求x=x2的解。

程序:#include "stdio.h"#include<math.h>main(){double x,xl,y,yl;int i,j;x=0.5;xl=x;y=0.5;yl=y;for(i=0;;i++){x=x*x;if(fabs(xl-x)<0.0001)break;else xl=x;}for(j=0;;j++){y=sqrt(y);if(fabs(yl-y)<0.0001)break;else yl=y;printf("x=%f,y=%f\n",x,y);}结果:x=0.000000,y=0.999915(五)牛顿迭代法y=f(x),求f(x*)=0。

数值分析(hilbert矩阵)病态线性方程组的求解matlab程序

数值分析(hilbert矩阵)病态线性方程组的求解matlab程序

(Hilbert 矩阵)病态线性方程组的求解理论分析表明,数值求解病态线性方程组很困难。

考虑求解如下的线性方程组的求解Hx = b ,期中H 是Hilbert 矩阵,()ij n n Hh ,11ij h i j ,i ,j = 1,2,…,n 1.估计矩阵的2条件数和阶数的关系2.对不同的n ,取(1,1,,1)nx K ?,分别用Gauss 消去,Jacobi 迭代,Gauss-seidel 迭代,SOR 迭代和共轭梯度法求解,比较结果。

3.结合计算结果,试讨论病态线性方程组的求解。

第1小题:condition.m %第1小题程序t1=20;%阶数n=20x1=1:t1;y1=1:t1;for i=1:t1H=hilb(i);y1(i)=log(cond(H));endplot(x1,y1);xlabel('阶数n');ylabel('2-条件数的对数(log(cond(H))');title('2-条件数的对数(log(cond(H))与阶数n 的关系图');t2=200;%阶数n=200x2=1:t2;y2=1:t2;for i=1:t2H=hilb(i);y2(i)=log(cond(H));endplot(x2,y2);xlabel('阶数n');ylabel('2-条件数的对数(log(cond(H))');title('2-条件数的对数(log(cond(H))与阶数n 的关系图');画出Hilbert 矩阵2-条件数的对数和阶数的关系n=200时n=20时从图中可以看出,1)在n小于等于13之前,图像近似直线log(cond(H))~1.519n-1.8332)在n大于13之后,图像趋于平缓,并在一定范围内上下波动,同时随着n的增加稍有上升的趋势第2小题:solve.m%m第2小题主程序N=4000;xGauss=zeros(N,1);xJacobi=zeros(N,1);xnJ=zeros(N,1);xGS=zeros(N,1);xnGS=zeros(N,1);xSOR=zeros(N,1);xnSOR=zeros(N,1);xCG=zeros(N,1);xnCG=zeros(N,1);for n=1:N;x=ones(n,1);t=1.1;%初始值偏差x0=t*x;%迭代初始值e=1.0e-8;%给定的误差A=hilb(n);b=A*x;max=100000000000;%可能最大的迭代次数w=0.5;%SOR迭代的松弛因子G=Gauss(A,b);[J,nJ]=Jacobi(A,b,x0,e,max);[GS,nGS]=G_S(A,b,x0,e,max);[S_R,nS_R]=SOR(A,b,x0,e,max,w);[C_G,nC_G]=CG(A,b,x0,e,max);normG=norm(G'-x);xGauss(n)=normG;normJ=norm(J-x);nJ;xJacobi(n)=normJ;xnJ(n)=nJ;normGS=norm(GS-x);nGS;xGS(n)=normGS;xnGS(n)=nGS;normS_R=norm(S_R-x);nS_R;xSOR(n)=normS_R;xnSOR(n)=nS_R;normC_G=norm(C_G-x);nC_G;xCG(n)=normC_G;xnCG(n)=nC_G;endGauss.m%Gauss消去法function x=Gauss(A,b)n=length(b);l=zeros(n,n);x=zeros(1,n);%消去过程for i=1:n-1for j=i+1:nl(j,i)=A(j,i)/A(i,i);for k=i:nA(j,k)=A(j,k)-l(j,i)*A(i,k);endb(j)=b(j)-l(j,i)*b(i);endend%回代过程x(n)=b(n)/A(n,n);for i=n-1:-1:1c=A(i,:).*x;x(i)=(b(i)-sum(c(i+1:n)))/A(i,i);endJacobi.m%Jacobi迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m 可能最大的迭代次数function [x,n]=Jacobi(A,b,x0,e,m)n=length(A);D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B=D\(L+U);f=D\b;x=B*x0+f;n=1;while norm(x-x0)>ex0=x;x=B*x0+f;n=n+1;if n>mdisp('Jacobi迭代次数过多,迭代可能不收敛');break;endendG_S.m%Gauss-Seidel迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数function [x,n]=G_S(A,b,x0,e,m)n=length(A);D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B=(D-L)\U;f=(D-L)\b;x=B*x0+f;n=1;while norm(x-x0)>ex0=x;x=B*x0+f;n=n+1;if n>mdisp('Gauss-Seidel迭代次数过多,迭代可能不收敛');break;endendSOR.m%SOR超松弛迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数,w松弛因子function [x,n]=SOR(A,b,x0,e,m,w)n=length(A);D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B=(D-w*L)\((1-w)*D+w*U);f=(D-w*L)\b*w;x=B*x0+f;n=1;while norm(x-x0)>ex0=x;x=B*x0+f;n=n+1;if n>mdisp('SOR超松弛迭代次数过多,迭代可能不收敛');break;endendCG.m%CG共轭梯度法,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数function [x,n]=CG(A,b,x0,e,m)r=b-A*x0;p=r;alpha=(r'*r)/(p'*(A*p));x=x0+alpha*p;r1=b-A*x;n=1;while norm(r1)>ebelta=(r1'*r1)/(r'*r);p=r1+belta*p;r=r1;x0=x;alpha=(r'*r)/(p'*(A*p));x=x0+alpha*p;r1=b-A*x;n=n+1;if n>mdisp('CG共轭梯度法迭代次数过多,迭代可能不收敛');break;endend。

数值分析中求解非线性方程的MATLAB求解程序

数值分析中求解非线性方程的MATLAB求解程序

数值分析中求解非线性方程的MATLAB求解程序1. fzero函数:fzero函数是MATLAB中最常用的求解非线性方程的函数之一、它使用了割线法、二分法和反复均值法等多种迭代算法来求解方程。

使用fzero函数可以很方便地求解单变量非线性方程和非线性方程组。

例如,要求解方程f(x) = 0,可以使用以下语法:``````2. fsolve函数:fsolve函数是MATLAB中求解多维非线性方程组的函数。

它是基于牛顿法的迭代算法来求解方程组。

使用fsolve函数可以非常方便地求解非线性方程组。

例如,要求解方程组F(x) = 0,可以使用以下语法:``````3. root函数:root函数是MATLAB中求解非线性方程组的函数之一、它采用牛顿法或拟牛顿法来求解方程组。

使用root函数可以非常方便地求解非线性方程组。

例如,要求解方程组F(x) = 0,可以使用以下语法:``````4. vpasolve函数:vpasolve函数是MATLAB中求解符号方程的函数。

它使用符号计算的方法来求解方程,可以得到精确的解。

vpasolve函数可以求解多变量非线性方程组和含有符号参数的非线性方程。

例如,要求解方程组F(x) = 0,可以使用以下语法:```x = vpasolve(F(x) == 0, x)```vpasolve函数会返回方程组的一个精确解x。

5. fsolve和lsqnonlin结合:在MATLAB中,可以将求解非线性方程转化为求解最小二乘问题的形式。

可以使用fsolve函数或lsqnonlin函数来求解最小二乘问题。

例如,要求解方程f(x) = 0,可以将其转化为最小二乘问题g(x) = min,然后使用fsolve或lsqnonlin函数来求解。

具体使用方法可以参考MATLAB官方文档。

6. Newton-Raphson法手动实现:除了使用MATLAB中的函数来求解非线性方程,还可以手动实现Newton-Raphson法来求解。

数值分析算法C语言程序

数值分析算法C语言程序

拉格朗日插值#in clude<stdio.h>#in clude<stdlib.h>#in clude<math.h>void Lagra nge(float s){double x[5]={0.2,0.4,0.6,0.8,1.0},y[5]={0.98,0.92,0.81,0.64,0.38},f,L=0; int i,j;for (i=0;i<5;i++){ f=1; for (j=0;j<5;j++) if(j!=i) f=(s-x[j])/(x[i]-x[j])*f; L+=f*y[i]; } printf("输出:%f\n",L);}void mai n(){float x;printf("输入插值点:");scan f("%f", &x);Lagra nge(x);}妙人斑值原:0-5爺出:0-871406any key to continue牛顿插值#in clude<stdlib.h>#in clude<stdio.h>#in clude<math.h> int ND(float s){double x[ 5]={0.2,0.4,0.6,0.8,1.0},y[5]={0.98,0.92,0.81,0.64,0.38},p=0,g,f;int i,j,k;for (i=0;i<5;i++){for (j=4;j>i;j--){ f=x[j]-x[j-i-1];y[j]=(y[j]-y[j-1])/f;}g=y[i+1];for (k=0;k<=i;k++)g=g*(s-x[k]);p=p+g;}printf("输出插值点函数值:%f\n",p+y[0]);return 1;void mai n(){float x;printf("输入插值点:"); scan f("%f", &x);ND(x);}三、埃尔米特插值#in clude<stdio.h>#in clude<stdlib.h>#in clude<math.h>void Hermite(float s){double x[3]={0.25,1,2.25},y[3]={0.125,1,3.375},z[3]={0.75,1.5,2.25};double H=0.0,a,b,f,g;int i,j;for (i=0;i<3;i++){f=1.0;g=0.0;for (j=0;j<3&&j!=i;j++) {f=f*(s-x[j])/(x[i]-x[j]);g=g+1/(x[i]-x[j]);} a=(1-2*(s-x[i])*g)*f*f;b=(s-x[i])*f*f;H=H+y[i]*a+z[i]*b;}prin tf("%f\n",H);}void mai n(){float x;printf("输入插值点:");scan f("%f", &x);Hermite(x);}四、三次样条插值#i nclude <math.h>#i nclude <stdio.h>#in clude <stdlib.h>void main(){int N=7,R=2,i,k;double p1,p2,p3,p4;double x[8]={0.5,0.7,0.9,1.1,1.3,1.5,1.7,1.9};double y[8]={0.4794,0.6442,0.7833,0.8912,0.9636,0.9975,0.9917,0.9463};double P0=-0.4794,Pn=-0.9463,u[3]={0.6,0.8,1.2},s[3];double h[7],a[8],c[7], g[8],af[8],ba[7],m[8];for(k=0;k <N;k++)h[k]=x[k+1]-x[k];for(k=1;k <N;k++) for(k=1;k <N;k++) for(k=1;k <N;k++) c[0]=a[N]=1; a[k]=h[k]/(h[k]+h[k-1]);c[k]=1-a[k];g[k]=3*(c[k]*(y[k+1]-y[k])/h[k]+a[k]*(y[k]-y[k-1])/h[k-1]);g[0]=3*(y[1]-y[0])/h[0]-P0*h[0]/2;g[N]=3*(y[N]-y[N-1])/h[N-1]+Pn*h[N-1]/2;ba[0]=c[0]/2; g[0]=g[0]/2; for(i=1;i <N;i++){ af[i]=2-a[i]*ba[i-1]; g[i]=(g[i]-a[i]*g[i-1])/af[i]; ba[i]=c[i]/af[i]; } af[N]=2-a[N]*ba[N-1]; g[N]=(g[N]-a[N]*g[N-1])/af[N]; m[N]=g[N];for(i=N-1;i>=0;i--) m[i]=g[i]-ba[i]*m[i+1];for(i=0;i <=R;i++){ k=0;while(u[i]> x[k+1])k++;p1=(h[k]+2*(u[i]-x[k])*pow((u[i]-x[k+1]),2)*y[k])/pow(h[k],3); p2=(h[k]-2*(u[i]-x[k+1])*pow((u[i]-x[k]),2)*y[k+1])/pow(h[k],3); p3=(u[i]-x[k])*pow((u[i]-x[k+1]),2)*m[k]/pow(h[k],2); p4=(u[i]-x[k+1])*pow((u[i]-x[k]),2)*m[k+1]/pow(h[k],2); s[i]=p1+p2+p3+p4;}printf( "\nx= ");for(i=0;i <=N;i++)printf( "%8.1f ",x[i]);printf( "\ny= ");for(i=0;i <=N;i++)printf( "%8.4f ",y[i]);printf( "\n\nu= ");for(i=0;i <=R;i++)printf( "%9.2f ",u[i]);printf( "\n 插值点:s= ");for(i=0;i <=R;i++)prin tf( "%9.5f ”,s[i]);prin tf("\n");五、复合梯形公式#in clude<stdio.h>#in clude<stdlib.h>#in clude<math.h> double FTX(i nt n,float a,float b){double f=0,t,h,*x,*y;int i;x=(double*)malloc(( n+1)*sizeof(double)); y=(double*)malloc(( n+1)*sizeof(double)); h=(b-a)/n;for(i=0;i< n+1;i++){x[i]=a+i*h;y[i]=s in (x[i]);}for(i=1;i< n; i++) f=f+2*y[i];t=h/2*(y[0]+f+y[ n]);printf("输出函数值:%f\n",t);return 1;}void mai n(){float a,b;int n;printf("输入区间上,下限:"); scan f("%f %f", &a, &b);printf("输入等分区间数:");scan f("%d",&n);FTX( n,a,b);}六、复合辛普森求积公式#in clude<stdio.h>#in clude<stdlib.h>#in clude<math.h> double FSP(int n,float a,float b){double f1=0,f2=0,h,*x1,*y1,*x2,*y2; int i;x1= (double*)malloc(( n+1)*sizeof(double));y1= (double*)malloc(( n+1)*sizeof(double));x2=(double*)malloc( n*sizeof(double));y2=(double*)malloc( n*sizeof(double)); h=(b-a)/n;for(i=0;i< n+1;i++) {x1[i]=a+h*i;y1[i]=si n(x1[i])*x1[i];} for(i=0;i< n;i++){x2[i]=x1[i]+h/2;y2[i]=si n( x2[i])*x2[i];}for(i=1;i< n;i++) f仁f1+2*y1[i];for(i=0;i< n;i++) f2=f2+4*y2[i];printf(” 输出函数值:%f\n",h/6*(y1[0]+f1+f2+y1[n]));return 1;}void mai n(){float a,b;int n;printf("输入区间上,下限:");scan f("%f %f", &a, &b);printf("输入等分区间数:");scan f("%d",&n);FSP( n,a,b);}七、直接三角分解法#in clude<stdio.h>#in clude<stdlib.h>#in clude<math.h>void mai n(){doubleA[3][3]={0.25,0.2,0.166667,0.3333,0.25,0.2,0.5,1,2},x[3],y[3],b[3]={9,8,8},L[3][3],U[3][3],f1=0, f2=0; int i,j,k;for(i=0;i<3;i++)for(j=0;j<3;j++){U[i][j]=O;L[i][j]=O;}for(i=0;i<3;i++){U[0][i]=A[0][i];L[i][0]=A[i][0]/U[0][0];L[i][i]=1;}for(i=1;i<3;i++)for(j=i;j<3;j++){for(k=0;k<=i-1;k++){f 1= f1+L[i][k]*U[k][j];f2+=L[j][k]*U[k][i];}U[i][j]=A[i][j]-f1;L[j][i]=(A[j][i]-f2)/U[i][i];f1=0;f2=0;}y[O]=b[O];for(i=1;i<3;i++){for(j=0;j<=i-1;j++) f1+=L[i][j]*y[j];y[i]=b[i]-f1;f1=0;}x[2]=y[2]/U[2][2];for(i=1;i>=0;i__){for(j=i+1;j<3;j++) f2+=U[i][j]*x[j];x[i]=(y[i]-f2)/U[i][i];f2=0;} printf(” 输出L 矩阵:\n”);for(i=0;i<3;i++){for(j=0;j<3;j++)prin tf("%f ",L[i][j]);prin tf("\n");}printf(” 输出U 矩阵:\n");for(i=0;i<3;i++){for(j=0;j<3;j++)prin tf("%f ",U[i][j]);prin tf("\n");}printf(”输出求解结果:\n”);for(i=0;i<3;i++)prin tf("%f ",x[i]);prin tf("\n");八、改进的平方法#in clude<stdio.h>#in clude<stdlib.h>#in clude<math.h>void mai n(){double A[3][3]={2,-1,1,-1,-2,3,1,3,1},x[3],y[3],b[3]={4,5,6},d[3 ],L [3][3],U[3][3],f1=0,f2=0;int i,j,k, n=3;for(i=0;i< n;i++)for(j=0;j< n;j++) {U[i][j]=0;L[i][j]=0;}d[0]=A[0][0];L[0][0]=1;for(i=1;i< n; i++){L[i][i]=1;for(j=0;j<=i-1;j++) {for(k=0;k<=j-1;k++) f1+=U[i][k]*L[j][k];U[i][j]=A[i][j]-f1;L[i][j]=U[i][j]/d[j];f2+=U[i][j]*L[i][j];f1=0;} d[i]=A[i][j]-f2;f2=0;} y[0]=b[0];for(i=1;i< n; i++) {for(j=0;j<=i-1;j++) f1+=L[i][j]*y[j];y[i]=b[i]-f1;f1=0;}x[ n-1]=y[ n-1]/d[ n-1];for(i=n-2;i>=0;i--) {for(j=i+1;j< n;j++) f2+=L[j][i]*x[j];x[i]=y[i]/d[i]-f2;f2=0;}printf(” 输出L 矩阵:\n");for(i=0;i< n;i++){for(j=0;j< n;j++) prin tf("%f ",L[i][j]); prin tf("\n");}printf("输出U 矩阵:\n");for(i=0;i< n;i++){for(j=0;j< n;j++) prin tf("%f ",U[i][j]); prin tf("\n");}printf("输出求解结果:\n");for(i=0;i< n;i++) prin tf("%f ",x[i]);prin tf("\n");}九、追赶法double A[5][5]={2,-1,0,0,0,-1,2,-1,0,0,0,-1,2,-1,0,0,0,-1,2,-1,0,0,0,-1,2},f[5]={1,0,0,0,0};double x[5],y[5],U[5][5];int i,j, n=5;for(i=0;i< n;i++)for(j=0;j<n;j++) U[i][j]=0;U[0][0]=U[ n-1][ n-1]=1;U[0][1]=A[0][1]/A[0][0];for(i=1;i< n-1;i++){ U[i][i]=1;U[i][i+1]=A[i][i+1]/(A[i][i]-A[i][i-1]*U[i-1][i]);}y[0]=f[0]/A[0][0];for(i=1;i< n; i++)y[i]=(f[i]-A[i][i-1]*y[i-1])/(A[i][i]-A[i][i-1]*U[i-1][i]);x[ n_1]=y[ n-1];for(i=n-2;i>=0;i--) x[i]=y[i]-U[i][i+1]*x[i+1];printf("输出U 矩阵:\n");for(i=0;i< n;i++){ for(j=0;j< n;j++) prin tf("%f ",U[i][j]); prin tf("\n");}printf(”输出求解结果:\n");for(i=0;i< n;i++) prin tf("%f ",x[i]);prin tf("\n");}十、雅可比迭代法#in clude<stdio.h>#in clude<stdlib.h>#in clude<math.h>void mai n(){double A[3][3]={5,2,1,-1,4,2,2,-3,10},x[50][3],b[3]={-12,20,3},f=0,L=1;int n=3,i,j,k=1;for (i=0;i<3;i++) x[0][i]=0;/* 初始值*/while(L<0.003){for(i=0;i <n ;i++){ for(j=0;j< n;j++)if(j!=i) f=f+A[i][j]*x[k-1][j];x[k][i]=(b[i]-f)/A[i][i];}L=x[k][0]-x[k-1][0];for(i=1;i <n ;i++)if((x[k][i]-x[k-1][i])/(x[k][i-1]-x[k-1][i-1])>1||(x[k][i]-x[k-1][i])/(x[k][i-1]-x[k-1][i-1])<-1) L=x[k][i]-x[k-1][i];k++;}printf(”输出每次迭代求得的X值:\n");for(i=1;i<k;i++){printf("第%d 次迭代:”,i);for(j=0;j< n;j++)prin tf("%f ",x[i][j]);prin tf("\n ”);}printf("\n 输出迭代次数:%d\n",k-1);咼斯一塞德尔迭代法double A[3][3]={5,2,1,-1,4,2,2,-3,10},x[50][3],b[3]={-12,20,3},f1=0,f2=0, L=1;int i,j,k=1, n=3;for (i=0;i<3;i++) x[0][i]=0;while(L>0.00001||L<-0.00001){for(i=0;i< n;i++){ for(j=0;j< n;j++){if(j<i) f1+=A[i][j]*x[k][j];else if(j>i) f2+=A[i][j]*x[k-1][j];}x[k][i]=(b[i]-f1-f2)/A[i][i];}L=x[k][0]-x[k-1][0];for(i=1;i <n ;i++)if((x[k][i]-x[k-1][i])/(x[k][i-1]-x[k-1][i-1])>1||(x[k][i]-x[k-1][i])/(x[k][i-1]-x[k-1][i-1])<-1) L=x[k][i]-x[k-1][i];k++;}printf(”输出迭代X矩阵:\n"); for(i=1;i<k;i++){for(j=0;j< n;j++) prin tf("%f ",x[i][j]);prin tf("\n");}printf("\n 输出迭代次数:%d\n",k-1);double A[3][3]={5,2,1,-1,4,2,2,-3,10},x[30][3],b[3]={-12,20,3},f1=0,f2=0, L=1,w=0.9; inti,j,k=1, n=3;for(i=0;i<3;i++) x[0][i]=0;while(L>0.0001||L<-0.0001){for(i=0;i <n ;i++){ for(j=0;j< n;j++){if(j<i) f1+=A[i][j]*x[k][j];else if(j>i) f2+=A[i][j]*x[k-1][j];}x[k][i]=w*(b[i]-f1-f2)/A[i][i];}L=x[k][0]-x[k-1][0];for(i=1;i< n; i++) if((x[k][i]-x[k-1][i])/(x[k][i-1]-x[k-1][i-1])>1||(x[k][i]-x[k-1][i])/(x[k][i-1]-x[k-1][i-1])<-1)L=x[k][i]-x[k-1][i];k++;}printf(”输出迭代X矩阵:\n");for(i=1;i<k;i++){for(j=0;j< n;j++) prin tf("%f ",x[i][j]); prin tf("\n");}printf("\n 输出迭代次数:%d\n",k-1);}数值分析算法程序班级:计算08Q2 班姓名:甄彦福。

数值分析 迭代法 C++程序

数值分析 迭代法 C++程序

课题三解线性方程组的迭代法实验目标:分别采用Jacobi迭代法,Gauss-Seidel迭代法和SOR迭代法求解线性方程组。

Jocabi迭代法:#include<iostream>#include<math.h>using namespace std;int i,j,k; //计数器int M = 2000;int Array(double ***Arr, int n){double **p;int i;p=(double **)malloc(n*sizeof(double *));if(!p)return 0;for(i=0;i<n;i++){p[i]=(double *)malloc(n*sizeof(double));if(!p[i])return 0;}*Arr=p;return 1;}void main(){double eps ;cout<<"默认最多迭代次数为2000次"<<endl<<"迭代精度为:";cin>>eps;double **matrix;int n;cout<<"矩阵大小为:";cin>>n;double *X;X= new double[n];double *Y;Y= new double[n];double *G;G= new double[n];for(i=0;i<n;i++){Y[i]=0;}if(!Array(&matrix,n))cout<<"内存分配失败!";elsecout<<"请输入矩阵:"<<endl;for( i=0;i<n;i++){for( j=0;j<n;j++){cin>>matrix[i][j];}}cout<<"请输入右端项:"<<endl;double *B;B = new double[n];for(i=0;i<n;i++){cin>>B[i];}for (i = 0 ;i< n;i++){if (fabs(matrix[i][i]) < eps){cout <<"打印失败"<<endl;return;}double T = matrix[i][i];for ( j = 0 ; j< n;j++){matrix[i][j] = -matrix[i][j]/T;}matrix[i][i] = 0;G[i] = B[i]/T;}int counter = 0;while (counter < M){for (i = 0;i < n; i++){double temp = 0;for (j = 0;j<n; j++){temp = temp + matrix[i][j]*Y[j];}X[i] = G[i] + temp;}double temp = 0;for (i = 0 ;i< n ; i++){temp = temp + fabs(X[i] - Y[i]);}if (temp <= eps)break;else{for( i = 0; i < n ;i++){Y[i] = X[i];}}counter++;}cout << "迭代次数为:"<<counter<<"次。

数值分析实验程序C语言

数值分析实验程序C语言
}
Y=Y+r*y[i];}
printf("%lf\n",Y);
}
2牛顿插值:
#include<stdio.h>
#define N 6
void main()
{
int i,k,n=N-1;
float X,Y,x[N],y[N],c[N],b[N];
printf("请输入xi:\n");
for(i=0;i<N;i++)
int k,N;
double x0,x1;
printf("请输入x0,N\n");
scanf("%lf,%d",&x0,&N);
for(k=1;k<=N;k++)
{
if(fd(x0)==0){printf("奇异标志\n");break;}
x1=x0-(f(x0)/fd(x0));
printf("X%d= %lf\n",k,x1);
}
}
float f(float x,float y)
{
float s;
s=y-(2*x)/y;
return s;
}
5牛顿迭代法:
#include<stdio.h>
#include<math.h>
#define e 1e-7
void main()
{
double f(double x);
double fd(double x);
1拉格朗日插值:
#include<stdio.h>
#define N 3

数值分析算法C语言程序

数值分析算法C语言程序

数值分析算法C语言程序数值分析是研究数学问题的近似解法的一门学科,其中包括了各种数值方法和算法。

本文将介绍数值分析中的常见算法,并给出相应的C语言程序。

1.二分法(Bisection Method)二分法是一种求函数零点的简单且常用的方法。

该方法的基本思想是通过不断将区间进行二分,并比较中点处函数值的正负来找到零点所在的区间。

```c#include <stdio.h>double f(double x)return x * x - 2;double bisection(double a, double b, double eps)double c;while ((b - a) > eps)c=(a+b)/2;if (f(c) == 0)break;}else if (f(a) * f(c) < 0)b=c;}elsea=c;}}return c;int maidouble a = 0.0;double b = 2.0;double result = bisection(a, b, eps);printf("The root is: %lf\n", result);return 0;```2.牛顿迭代法(Newton's Method)牛顿迭代法是一种高效的求函数零点的方法。

该方法的基本思想是通过对函数进行线性逼近,不断逼近函数的零点。

```c#include <stdio.h>#include <math.h>double f(double x)return x * x - 2;double df(double x)return 2 * x;double newton(double x0, double eps) double x = x0;double deltaX = f(x) / df(x);while (fabs(deltaX) > eps)deltaX = f(x) / df(x);x = x - deltaX;}return x;int maidouble x0 = 2.0;double result = newton(x0, eps); printf("The root is: %lf\n", result); return 0;```3.高斯消元法(Gaussian Elimination)高斯消元法是一种用于求解线性方程组的方法。

数值分析中求解非线性方程的MATLAB求解程序(6种)

数值分析中求解非线性方程的MATLAB求解程序(6种)

数值分析中求解非线性方程的MATLAB求解程序(6种)数值分析中求解非线性方程的MATLAB求解程序(6种)1.求解不动点function [k,p,err,P]=fixpt(g,p0,tol,max1)%求解方程x=g(x) 的近似值,初始值为p0%迭代式为Pn+1=g(Pn)%迭代条件为:在迭代范围内满足|k|<1(根及附近且包含初值)k为斜率P(1)=p0;for k=2:max1P(k)=feval(g,P(k-1));err=abs(P(k)-P(k-1));relerr=err/(abs(P(k))+eps);p=P(k);if (err<tol)|(relerr<tol)< p="">break;endendif k==max1disp('超过了最长的迭代次数')endP=P';2.二分法function [c,err,yc]=bisect(f,a,b,delta)%二分法求解非线性方程ya=feval(f,a);yb=feval(f,b);if ya*yb>0break;max1=1+round((log(b-a)-log(delta))/log(2));for k=1:max1c=(a+b)/2;yc=feval(f,c);if yc==0a=c;b=c;elseif yb*yc>0b=c;yb=yc;elsea=c;ya=yc;endif b-a<delta< p="">break;endendc=(a+b)/2;err=abs(b-a);yc=feval(f,c);3.试值法function [c,err,yc]=regula(f,a,b,delta,epsilon,max1) %试值法求解非线性方程%f(a)和飞(b)异号ya=feval(f,a);yb=feval(f,b);if ya*yb>0disp('Note:f(a)*f(b)>0');for k=1:max1dx=yb*(b-a)/(yb-ya);c=b-dx;ac=c-a;yc=feval(f,c);if yc==0break;elseif yb*yc>0b=c;yb=yc;elsea=c;ya=yc;enddx=min(abs(dx),ac);if abs(dx)<delta|abs(yc)<epsilon< p="">break;endendc;err=abs(b-a)/2;yc=feval(f,c);4.求解非线性方程根的近似位置function R=approot(X,epsilon)%求解根近似位置%为了粗估算方程f(x)=0在区间[a,b]的根的位置,%使用等间隔采样点(xk,f(xk))和如下的评定准则:%f(xk-1)与f(xk)符号相反,%或者|f(xk)|足够小且曲线y=f(x)的斜率在%(xk,f(xk))附近改变符号。

数值分析实验,用程序实现Hermite插值法

数值分析实验,用程序实现Hermite插值法

《数值分析》实验报告实验序号:实验六 实验名称: Hermite 插值法1. 实验目的:学会Hermite 插值法,并应用该算法于实际问题.2. 实验内容:求一个函数ϕ(x )用来近似函数f (x ),用分段三次Hermit 插值的方法来求解近似函数ϕ(x )并画出近似函数图像及原函数图像。

设在区间[a,b]上,给定n+1个插值节点b x x x x a n =<<<<=...210和相应的函数值n y y y ,...,,10以及一阶导数值''1'0,...,,ny y y ,求一个插值函数)(x H ,满足以下条件: (1)),...,2,1,0()(,)(''n i y x H y x H i i i i === (2) )(x H 在每一个小区间[1,+j j x x ]上是三次多项式。

对于给定函数11-,2511)(2≤≤+=x x x f 。

在区间[]11-,上画出f (x )和分段三次Hermit 插值函数)(x H 的函数图像。

3. 实验分析:算法分析:1. 分段三次Hermit 插值的算法思想:分段三次Hermit 插值的做法是在每一个小区间上作三次Hermit 插值,因此在每一个插值节点上都需要构造两个插值基函数)(),(x H x h i i ,然后再作它们的线性组合。

分段三次Hermit 插值基函数如下:⎪⎩⎪⎨⎧≤≤----+=其它 0 ))(21()(1021010100x x x x x x x x x x x x h ⎪⎩⎪⎨⎧≤≤---=其它 0 ))(()(10210100x x x x x x x x x x H1,...,2,1 0 ))(21( ))(21()(1211112111-=⎪⎪⎪⎩⎪⎪⎪⎨⎧≤<----+≤≤----+=++++---n i x x x x x x x x x x x x x x x x x x x x x x x h i i i i i i i i i i-i i i i i i i 其它1,...,2,1 0 ))(( ))(()(12111211-=⎪⎪⎪⎩⎪⎪⎪⎨⎧≤<---≤≤---=+++--n i x x x x x x x x x x x x x x x x x x x H i i i i i i i i-i i i i i 其它⎪⎩⎪⎨⎧≤<----+=---其它 0 ))(21()(1-n 2111n n n n n n n n x x x x x x x x x x x x h ⎪⎩⎪⎨⎧≤<---=--其它 0 ))(()(1-n 211n n n n n n x x x x x x x x x x H 分段三次Hermit 插值函数是:∑=+=n i i i i i x H y x h y x H 0'))()(()( 4. 实验代码:// LDlg.cpp : implementation file//#include "stdafx.h"#include "L.h"#include "LDlg.h"#ifdef _DEBUG#define new DEBUG_NEW#undef THIS_FILEstatic char THIS_FILE[] = __FILE__;#endif/////////////////////////////////////////////////////////////////////////////// CAboutDlg dialog used for App Aboutclass CAboutDlg : public CDialog{public:CAboutDlg();// Dialog Data//{{AFX_DATA(CAboutDlg)enum { IDD = IDD_ABOUTBOX };//}}AFX_DATA// ClassWizard generated virtual function overrides//{{AFX_VIRTUAL(CAboutDlg)protected:virtual void DoDataExchange(CDataExchange* pDX); // DDX/DDV support //}}AFX_VIRTUAL// Implementationprotected://{{AFX_MSG(CAboutDlg)//}}AFX_MSGDECLARE_MESSAGE_MAP()};CAboutDlg::CAboutDlg() : CDialog(CAboutDlg::IDD){//{{AFX_DATA_INIT(CAboutDlg)//}}AFX_DATA_INIT}void CAboutDlg::DoDataExchange(CDataExchange* pDX){CDialog::DoDataExchange(pDX);//{{AFX_DATA_MAP(CAboutDlg)//}}AFX_DATA_MAP}BEGIN_MESSAGE_MAP(CAboutDlg, CDialog)//{{AFX_MSG_MAP(CAboutDlg)// No message handlers//}}AFX_MSG_MAPEND_MESSAGE_MAP()///////////////////////////////////////////////////////////////////////////// // CLDlg dialogCLDlg::CLDlg(CWnd* pParent /*=NULL*/): CDialog(CLDlg::IDD, pParent){//{{AFX_DATA_INIT(CLDlg)// NOTE: the ClassWizard will add member initialization here//}}AFX_DATA_INIT// Note that LoadIcon does not require a subsequent DestroyIcon in Win32m_hIcon = AfxGetApp()->LoadIcon(IDR_MAINFRAME);}void CLDlg::DoDataExchange(CDataExchange* pDX){CDialog::DoDataExchange(pDX);//{{AFX_DATA_MAP(CLDlg)// NOTE: the ClassWizard will add DDX and DDV calls here//}}AFX_DATA_MAP}BEGIN_MESSAGE_MAP(CLDlg, CDialog)//{{AFX_MSG_MAP(CLDlg)ON_WM_SYSCOMMAND()ON_WM_PAINT()ON_WM_QUERYDRAGICON()ON_BN_CLICKED(IDC_LARGRI, OnLargri)ON_BN_CLICKED(IDC_BUTTON2, OnButton2)ON_BN_CLICKED(IDC_HERMITE, OnHermite)//}}AFX_MSG_MAPEND_MESSAGE_MAP()///////////////////////////////////////////////////////////////////////////// // CLDlg message handlersBOOL CLDlg::OnInitDialog(){CDialog::OnInitDialog();// Add "About..." menu item to system menu.// IDM_ABOUTBOX must be in the system command range.ASSERT((IDM_ABOUTBOX & 0xFFF0) == IDM_ABOUTBOX);ASSERT(IDM_ABOUTBOX < 0xF000);CMenu* pSysMenu = GetSystemMenu(FALSE);if (pSysMenu != NULL){CString strAboutMenu;strAboutMenu.LoadString(IDS_ABOUTBOX);if (!strAboutMenu.IsEmpty()){pSysMenu->AppendMenu(MF_SEPARATOR);pSysMenu->AppendMenu(MF_STRING, IDM_ABOUTBOX, strAboutMenu);}}// Set the icon for this dialog. The framework does this automatically // when the application's main window is not a dialogSetIcon(m_hIcon, TRUE); // Set big iconSetIcon(m_hIcon, FALSE); // Set small icon// TODO: Add extra initialization herereturn TRUE; // return TRUE unless you set the focus to a control}void CLDlg::OnSysCommand(UINT nID, LPARAM lParam){if ((nID & 0xFFF0) == IDM_ABOUTBOX){CAboutDlg dlgAbout;dlgAbout.DoModal();}else{CDialog::OnSysCommand(nID, lParam);}// If you add a minimize button to your dialog, you will need the code below // to draw the icon. For MFC applications using the document/view model, // this is automatically done for you by the framework.void CLDlg::OnPaint(){if (IsIconic()){CPaintDC dc(this); // device context for paintingSendMessage(WM_ICONERASEBKGND, (WPARAM) dc.GetSafeHdc(), 0);// Center icon in client rectangleint cxIcon = GetSystemMetrics(SM_CXICON);int cyIcon = GetSystemMetrics(SM_CYICON);CRect rect;GetClientRect(&rect);int x = (rect.Width() - cxIcon + 1) / 2;int y = (rect.Height() - cyIcon + 1) / 2;// Draw the icondc.DrawIcon(x, y, m_hIcon);}else{CDialog::OnPaint();}}// The system calls this to obtain the cursor to display while the user drags // the minimized window.HCURSOR CLDlg::OnQueryDragIcon(){return (HCURSOR) m_hIcon;}void CLDlg::OnOK(){int x00=300,y00=350,i,j;double x;CDC *pDC=GetDC();pDC->SetMapMode(MM_LOMETRIC);pDC->SetViewportOrg(x00,y00);//画坐标轴与原函数for(i=-700; i<=700; i++){pDC->SetPixel(i,0,RGB(0,0,0));pDC->SetPixel(0,i,RGB(0,0,0));}for(x=-1; x<=1; x+=0.001){double j=400.0/(1+25*x*x);pDC->SetPixel(x*500,j,RGB(255,0,0));}pDC->TextOut(-30,-10,"0");pDC->TextOut(-30,430,"1");pDC->TextOut(490,-10,"1");pDC->TextOut(-490,-10,"-1");pDC->MoveTo(-10,680); //x箭头pDC->LineTo(0,700);pDC->MoveTo(0,700);pDC->LineTo(10,680);pDC->MoveTo(680,10); //y箭头pDC->LineTo(700,0);pDC->MoveTo(700,0);pDC->LineTo(680,-10);pDC->TextOut(-30,700,"y");pDC->TextOut(700,-10,"x");}void CLDlg::OnLargri(){int x00=300,y00=350,i,j;double x;CDC *pDC=GetDC();pDC->SetMapMode(MM_LOMETRIC);pDC->SetViewportOrg(x00,y00);//画坐标轴for(i=-700; i<=700; i++){pDC->SetPixel(i,0,RGB(0,0,0));pDC->SetPixel(0,i,RGB(0,0,0));}double yx[]={-1,-0.8,-0.6,-0.4,-0.2,0,0.2,0.4,0.6,0.8,1}; pDC->TextOut(-30,-10,"0");pDC->TextOut(-30,430,"1");pDC->TextOut(490,-10,"1");pDC->TextOut(-490,-10,"-1");pDC->MoveTo(-10,680); //x箭头pDC->LineTo(0,700);pDC->MoveTo(0,700);pDC->LineTo(10,680);pDC->MoveTo(680,10); //y箭头pDC->LineTo(700,0);pDC->MoveTo(700,0);pDC->LineTo(680,-10);pDC->TextOut(-30,700,"y");pDC->TextOut(700,-10,"x");// 拉格朗日差值的函数double yy[12],lx[12],ly[12];double l_fenzi[12],l_fenmu[12];double l_x,l_y;for(i=0; i<=10; i++){yy[i]=1.0/(1+25*yx[i]*yx[i]);for(i=0; i<=10; i++){l_fenmu[i]=1.0;for(j=0; j<=10; j++){if(i!=j)l_fenmu[i]=l_fenmu[i]*(yx[i]-yx[j]);}}double qq,pp;for(qq=-1; qq<=1; qq+=0.0003){for(i=0; i<=10; i++){l_fenzi[i]=1.0;for(j=0; j<=10; j++){if(i!=j)l_fenzi[i]=l_fenzi[i]*(qq-yx[j]);}}pp=0;for(i=0; i<=11; i++){pp=pp+1.0/(1+25*yx[i]*yx[i])*l_fenzi[i]/l_fenmu[i];}pDC->SetPixel(qq*500,pp*390+5,RGB(132,112,225));}void CLDlg::OnButton2(){int x00=300,y00=350,i,j;double x;CDC *pDC=GetDC();pDC->SetMapMode(MM_LOMETRIC);pDC->SetViewportOrg(x00,y00);//画坐标轴与原函数for(i=-700; i<=700; i++){pDC->SetPixel(i,0,RGB(0,0,0));pDC->SetPixel(0,i,RGB(0,0,0));}double yx[]={-1,-0.8,-0.6,-0.4,-0.2,0,0.2,0.4,0.6,0.8,1}; double yy[14];for(i=0; i<=10; i++){yy[i]=1.0/(1+25*yx[i]*yx[i]);}pDC->TextOut(-30,-10,"0");pDC->TextOut(-30,430,"1");pDC->TextOut(490,-10,"1");pDC->TextOut(-490,-10,"-1");pDC->MoveTo(-10,680); //x箭头pDC->LineTo(0,700);pDC->MoveTo(0,700);pDC->LineTo(10,680);pDC->MoveTo(680,10); //y箭头pDC->LineTo(700,0);pDC->MoveTo(700,0);pDC->LineTo(680,-10);pDC->TextOut(-30,700,"y");pDC->TextOut(700,-10,"x");// 线性分段差值的图像CPen pen;CPen*oldpen;pen.CreatePen(PS_SOLID,5,RGB(0,0,0));oldpen=pDC->SelectObject(&pen);for(i=0; i<10; i++){pDC->MoveTo(yx[i]*480,yy[i]*400);pDC->LineTo(yx[i+1]*480,yy[i+1]*400); }}void CLDlg::OnHermite(){int x00=300,y00=350,i,j;double x;CDC *pDC=GetDC();pDC->SetMapMode(MM_LOMETRIC);pDC->SetViewportOrg(x00,y00);//画坐标轴与原函数for(i=-700; i<=700; i++){pDC->SetPixel(i,0,RGB(0,0,0));pDC->SetPixel(0,i,RGB(0,0,0));}double yx[]={-1,-0.8,-0.6,-0.4,-0.2,0,0.2,0.4,0.6,0.8,1};double yy[12];for(i=0; i<=10; i++){yy[i]=1.0/(1+25*yx[i]*yx[i]);}pDC->TextOut(-30,-10,"0");pDC->TextOut(-30,430,"1");pDC->TextOut(490,-10,"1");pDC->TextOut(-490,-10,"-1");pDC->MoveTo(-10,680); //x箭头pDC->LineTo(0,700);pDC->MoveTo(0,700);pDC->LineTo(10,680);pDC->MoveTo(680,10); //y箭头pDC->LineTo(700,0);pDC->MoveTo(700,0);pDC->LineTo(680,-10);pDC->TextOut(-30,700,"y");pDC->TextOut(700,-10,"x");//分段三次Hermite差值的函数double x0,x1,yd1,yd0,y1,y0;for(i=0; i<10; i++){x0=yx[i],x1=yx[i+1];y0=1.0/(1+25*x0*x0);y1=1.0/(1+25*x1*x1);yd0=-(50*x0)*1.0/((1+25*x0*x0)*(1+25*x0*x0));yd1=-(50*x1)*1.0/((1+25*x1*x1)*(1+25*x1*x1));for(double qq=x0; qq<x1; qq+=0.005){double pp= y0*(1+2*(qq-x0)/(x1-x0)) * (qq-x1)/(x0-x1) * (qq-x1)/(x0-x1)+y1*(1+2*(qq-x1)/(x0-x1)) * (qq-x0)/(x1-x0) * (qq-x0)/(x1-x0)+yd0*(qq-x0) * (qq-x1)/(x0-x1) * (qq-x1)/(x0-x1)+yd1*(qq-x1) * (qq-x0)/(x1-x0) * (qq-x0)/(x1-x0);pDC->SetPixel(qq*500,pp*400,RGB(225,185,15));}}}5.实验截图6. 实验结果分析:通过本次实验我对分段三次Hermit插值有了更深刻更全面的掌握,它在给定了节点处的函数值和导数值以后,构造了一个整体上具有一阶连续微商的插值函数。

数值分析(计算方法)介绍

数值分析(计算方法)介绍

Zeno悖论所描述的逼近过程正是这种迭代过程,当k→∞时,tk →t* (问题2: 证明该结论!)。大家知道,任何形式的重复都可看成是 “时间”的量度。Zeno在刻画人龟追赶问题中设置了两个“时钟”:一 个是日常的钟,另外Zeno又将迭代次数视为另一种时钟,不妨称之为 Zeno钟。Zeno公式(2)表明,当Zeno钟趋于∞时人才能追上龟,Zeno 正是据此断言人永远追不上龟。
9
North China Elec. P.U.
Numerical Analysis
2019/12/20
算法设计技术
J. G. Liu
引例
古希腊哲学家Zeno在两千多年前提出过一个骇人听闻的命题: 一个人不管跑得多快,也追不上爬在他前面的一只乌龟。这就 是著名的Zeno悖论。
Zeno在论证这个命题时采取了如下形式的逻辑推理:设人与龟 同时同向起跑,如果龟不动,那么人经过某段时间便能追上它; 但实际上在这段时间内龟又爬了一段路程,从而人又得重新追 赶,如下图所示,这样每追赶一次所归结的是同样类型的追赶 问题,因而这种追赶过程“永远”不会终结。
Numerical Analysis
2019/12/20
J. G. Liu
数值分析
——插值、拟合与数值微积分
主讲: 刘敬刚
School of Math. & Phys.
1
North China Elec. P.U.
Numerical Analysis
2019/12/20
数值分析(计算方法)简介
• 引例
School of Math. & Phys.
11
North China Elec. P.U.
Numerical Analysis

三种迭代法matlab程序 数值分析

三种迭代法matlab程序 数值分析


end
• end
• err=abs(norm(X'-P));
• P=X';
• if(err<delta)

break
• end
• end
• X=X';
• err,k
function X=jacobi(A,b,P,delta,max1) %A是n维非奇异阵。%b是n维向量。%P是初值。%delta是误差界。 %max1是给定的迭代最高次数。%X为所求的方程组AX=b的近似解。 N=length(b); for k=1:max1 for j=1:N
X(j)=(b(j)-A(j,[1:j-1,j+1:N])*P([1:j-1,j+1:N]))/A(j,j); end err=abs(norm(X'-P)); P=X'; if(err<delta)
break end end X=X';k,erction X=gseid(A,b,P,delta,max1)
• %A是n维非奇异阵。%b是n维向量。%P是初值。
• %delta是误差界。%max1是给定的迭代最高次数。%X为所求的方程组AX=b的近似解。
• N=length(b);
• for k=1:max1
• for j=1:N

if j==1

X(1)=(b(1)-A(1,2:N)*P(2:N))/A(1,1);

elseif j==N

X(N)=(b(N)-A(N,1:N-1)*(X(1:N-1))')/A(N,N);

else

X(j)=(b(j)-A(j,1:j-1)*X(1:j-1)-A(j,j+1:N)*P(j+1:N))/A(j,j);

数值分析各算法流程图

数值分析各算法流程图

01,,n1,,n1,,)n x及数值分析各算法流程图一、插值1、 拉格朗日插值流程图:( 相应程序:lagrintp(x,y,xx))2,,n ,,j n 1,2,,n 1,,)n 2、 牛顿插值流程图(1)产生差商表的算法流程图(相应程序:divdiff(x,y))注:1、另一程序divdiff1(x,y),输出的矩阵包含了节点向量。

而divdiff(x,y)不含节点向量。

2、另一程序tableofdd(x,y,m),输出的是表格形式,添加了表头。

1,,),,n m 及1,,m (2)非等距节点的牛顿插值流程图(相应程序:newtint11(x,y,xx,m)) 、注:1、虽然程序newtint11(x,y,xx,m)考虑了多种情形,看上去很复杂,但基本流程结构还是如上图所示。

2、程序中调用的子程序是divdiff 。

若调用的子程序是divdiff1的话,流程图中的第三,第四,第五步要相应的改一下数字。

2,3,,1m +1,,j1,2,,n=1,2,,)n m 及(3)求差分表的流程图(相应程序:difference(y,m))注:1、difference 输出的是矩阵D 。

而另一程序tableofd(y,m),输出的是带有表头的差分表。

n x m1,,),,1,,m注:1、程序newtforward1(x,y,xx,m))的结构与上述流程图一致,xx可以是数组。

2、另一程序newtforward(x,y,xx,m))先求出插值多项式,再求插值多项式在插值点的函数值。

基本结构还是和上面的流程图一样。

n x m1,,),,-x x1,,m注:1、程序newtbackward1(x,y,xx,m))的结构与上述流程图一致,xx可以是数组。

2、另一程序newtbackward(x,y,xx,m))先求出插值多项式,再求插值多项式在插值点的函数值。

基本结构还是和上面的流程图一样。

1,2,,n1,2,,n ,2,,)n x及3、Hermite 插值流程图(1) 已知条件中一阶导数的个数与插值节点的个数相等时的Hermite 插值流程图。

数值分析—二分法解非线性方程组—FORTRAN程序

数值分析—二分法解非线性方程组—FORTRAN程序

数值分析—二分法解非线性方程组—FORTRAN程序以下为一个使用二分法解非线性方程组的Fortran程序:```fortranPROGRAMBISECTION_METHODIMPLICITNONEINTEGER :: n, i, max_iterREAL :: a, b, tol, xn, fa, fb, fn!假设我们要解的方程组是f(x)=0,可以在这里定义f(x) REALFUNCTIONf(x)REAL::x!计算f(x)的值,并返回f=x**2-2ENDFUNCTIONf!主程序开始WRITE(*,*)'输入方程的左右端点a和b:'READ(*,*)a,bWRITE(*,*)'输入迭代的最大次数和误差容限:'READ(*,*) max_iter, tol!输出表头WRITE(*, '(3X, "Iteration", 3X, "a", 6X, "b", 6X, "f(a)", 4X, "f(b)", 4X, "f(n)")')WRITE(*,'(A)')'------------------------------------------------------------------'!二分法迭代DO i = 1, max_iter! 计算中点 xnxn = (a + b) / 2!计算f(a)、f(b)和f(n)fa = f(a)fb = f(b)fn = f(xn)!输出迭代步骤结果WRITE(*,'(I5,3X,F10.6,3X,F10.6,3X,F10.6,3X,F10.6,3X,F10.6)') & i, a, b, fa, fb, fn!判断迭代是否结束IF (ABS(fn) < tol) EXIT!根据f(a)、f(b)和f(n)判断下一步迭代的区间IF (fa * fn < 0) THENb = xnELSEa = xnENDIFENDDOWRITE(*,'(A)')'------------------------------------------------------------------'WRITE(*,*)'迭代结束!最终结果:'WRITE(*, '(A)') 'x =', xnENDPROGRAMBISECTION_METHOD```这个程序使用二分法迭代解决非线性方程组。

数值分析方法

数值分析方法

数值分析方法数值分析方法是一种利用数值计算来解决数学问题的方法。

它主要应用于工程、科学和其他领域中的实际问题,例如求解微分方程、插值、逼近、线性代数问题等。

数值分析方法的主要目标是通过数值计算来获得数学问题的近似解,因为很多实际问题并没有精确的解析解。

在本文中,我们将介绍数值分析方法的基本原理和常用技术,以及其在实际问题中的应用。

数值分析方法的基本原理是将连续的数学问题转化为离散的数值计算问题。

在实际应用中,我们通常会遇到无法通过解析方法求解的复杂问题,这时就需要借助数值分析方法来进行求解。

数值分析方法的主要步骤包括建立数学模型、离散化、选择适当的数值计算方法、计算和分析结果。

在建立数学模型时,我们需要将实际问题抽象为数学问题,并选择合适的数学方程描述。

离散化是将连续的数学问题转化为离散的数值计算问题,通常包括网格化、时间步长等处理。

选择适当的数值计算方法是数值分析方法中非常重要的一步,常用的数值计算方法包括有限差分法、有限元法、数值积分法等。

最后,我们需要进行数值计算,并对结果进行分析和验证。

数值分析方法在实际问题中有着广泛的应用。

在工程领域,数值分析方法被广泛应用于结构分析、流体力学、热传导等问题的求解。

在科学领域,数值分析方法被用于模拟天体运动、地球物理问题、量子力学等领域。

另外,在金融领域,数值分析方法也被用于期权定价、风险管理等问题的求解。

可以说,数值分析方法已经成为现代科学技术发展中不可或缺的工具之一。

总之,数值分析方法是一种通过数值计算来解决数学问题的方法,它的基本原理是将连续的数学问题转化为离散的数值计算问题。

在实际应用中,数值分析方法有着广泛的应用,包括工程、科学、金融等领域。

通过数值分析方法,我们可以求解很多无法通过解析方法求解的复杂问题,为科学技术的发展提供了重要的支持。

希望本文能够对读者对数值分析方法有所了解,并在实际问题中能够灵活运用数值分析方法来解决问题。

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

用N-C积分公式计算sin(x)在区间[0,∏]上的积分值。

#include"stdio.h"
#include"math.h"
void main()
{
int n,k;
double sum=0.0,a=0.0,b=3.2415926;
double Cotes[8][9]={{0.5,0.5},{1.0/6,4.0/6,1.0/6},{1.0/8,3.0/8,3.0/8,1.0/8},
{7.0/90,32.0/90,12.0/90,32.0/90,7.0/90},{19.0/288,75.0/288,50.0/288,50.0/288,75.0/288, 19.0/288},
{41.0/840,216.0/840,27.0/840,272.0/840,27.0/840,216.0/840,41.0/840},
{751.0/17280,3577.0/17280,1323.0/17280,2989.0/17280,2989.0/17280,1323.0/17280,35 77.0/17280,751.0/17280},
{989.0/28350,5888.0/28350,-928.0/28350,10496.0/28350,-4540.0/28350,10496.0/28350, -928.0/28350,5888.0/28350,989.0/28350}};
//printf("请输入积分区间a和b:");
//scanf("%lf,%lf",&a,&b);
printf("请输入积分节点n(1<=n<=8):");
scanf("%d",&n);
printf("\n");
for(k=0;k<=n;k++)
sum+=Cotes[n-1][k]*(sin(a+k*(b-a)/n));
sum=sum*(b-a);
printf("%lf\n",sum);
printf("误差值为:%lf\n",2.0-sum);
}
用复化梯形公式和复化Sipsion公式求函数(sinx)/x在区间[0,1]上的积分。

#include"stdio.h"
#include"math.h"
void Ti(double a,double b,int n)
{
double sum=0.0,h=0.0,x=0.0;
int i;
h=(b-a)/n;
for(i=1;i<n;i++)
{
x=a+i*h;
sum+=2*sin(x)/(x);
}
sum=(h/2)*(1+sin(b)/b+sum);
printf("函数sin(x)/x在区间[%lf,%lf]上的积分是%lf\n",a,b,sum);
}
void Si(double a,double b,int n)
{
double sum=0.0,h=0.0,x=0.0,y=0.0;
int i,j;
h=(b-a)/n;
for(i=1;i<=n/2;i++)
{
x=a+(2*i-1)*h;
sum+=4*sin(x)/x;
}
for(j=1;j<=n/2-1;j++)
{
y=a+2*j*h;
sum+=2*sin(y)/y;
}
sum=(h/3)*(1+sin(b)/b+sum);
printf("函数sin(x)/x在区间[%lf,%lf]上的积分是%lf\n",a,b,sum);
}
void main()
{
int s,n;
double a=0.0,b=0.0;
printf("请选择公式(1为复化梯形公式,2为复化Simpsion公式):");
scanf("%d",&s);
printf("\n");
printf("请输入分段数目n:");
scanf("%d",&n);
printf("\n");
printf("请输入求积区间a、b:");
scanf("%lf,%lf",&a,&b);
if(s==1)
Ti(a,b,n);//调用复化梯形公式else if(s==2)
Si(a,b,n);//调用复化Simpion公式}。

相关文档
最新文档