数值计算方法大作业

合集下载

数值分析大作业一

数值分析大作业一

数值分析大作业一一、算法设计方案1、求λ1和λ501的值:思路:采用幂法求出按模最大特征值λmax,该值必为λ1或λ501,若λmax小于0,则λmax=λ1;否则λmax=λ501。

再经过原点平移,使用幂法迭代出矩阵A-λmax I的特征值,此时求出的按模最大特征值即为λ1和λ501的另一个值。

2、求λs的值:采用反幂法求出按模最小的特征值λmin即为λs,其中的方程组采用LU分解法进行求解。

3、求与μk最接近的特征值:对矩阵A采用带原点平移的反幂法求解最小特征值,其中平移量为:μk。

4、A的条件数cond(A)=| λmax/λmin|;5、A的行列式的值:先将A进行LU分解,再求U矩阵对角元素的乘积即为A 行列式的值。

二、源程序#include<iostream>#include<iomanip>#include<math.h>#define N 501#define E 1.0e-12 //定义精度常量#define r 2#define s 2using namespace std;double a[N];double cc[5][N];void init();double mifa();double fmifa();int max(int aa,int bb);int min(int aa,int bb);int max_3(int aa,int bb,int cc);void LU();void main(){double a1,a2,d1,d501=0,ds,det=1,miu[39],lamta,cond;int i,k;init();/*************求λ1和λ501********************/a1=mifa();if(a1<0)d1=a1; //若小于0则表示λ1的值elsed501=a1; //若大于0则表示λ501的值for(i=0;i<N;i++)a[i]=a[i]-a1;a2=mifa()+a1;if(a2<0)d1=a2; //若小于0则表示λ1的值elsed501=a2; //若大于0则表示λ501的值cout<<"λ1="<<setiosflags(ios::scientific)<<setprecision(12)<<d1<<"\t";cout<<"λ501="<<setiosflags(ios::scientific)<<setprecision(12)<<d501<<endl;/**************求λs*****************/init();ds=fmifa();cout<<"λs="<<setiosflags(ios::scientific)<<setprecision(12)<<ds<<endl;/**************求与μk最接近的特征值λik**************/cout<<"与μk最接近的特征值λik:"<<endl;for(k=0;k<39;k++){miu[k]=d1+(k+1)*(d501-d1)/40;init();for(i=0;i<N;i++)a[i]=a[i]-miu[k];lamta=fmifa()+miu[k];cout<<"λi"<<k+1<<"\t\t"<<setiosflags(ios::scientific)<<setprecision(12)<<lamta<<en dl;}/**************求A的条件数**************/cout<<"矩阵A的条件式";cond=abs(max(abs(d1),abs(d501))/ds);cout<<"cond="<<setiosflags(ios::scientific)<<setprecision(12)<<cond<<endl;/**************求A的行列式**************/cout<<"矩阵A的行列式";init();LU();for(i=0;i<N;i++){det*=cc[2][i];}cout<<"det="<<setiosflags(ios::scientific)<<setprecision(12)<<det<<endl;system("pause");}/**************初始化函数,给a[N]赋值*************/void init(){int i;for(i=1;i<=501;i++)a[i-1]=(1.64-0.024*i)*sin((double)(0.2*i))-0.64*exp((double)(0.1/i)); }/**************幂法求最大绝对特征值**************/double mifa(){int i,k=0;double u[N],y[N]={0},b=0.16,c=-0.064,Beta_=0,error;for(i=0;i<501;i++)u[i]=1; //令u[N]=1for(k=1;k<2000;k++) //控制最大迭代次数为2000{/***求y(k-1)***/double sum_u=0,gh_sum_u;for(i=0;i<N;i++){sum_u+=u[i]*u[i]; }gh_sum_u=sqrt(sum_u);for(i=0;i<N;i++){y[i]=u[i]/gh_sum_u;}/****求新的uk****/u[0]=a[0]*y[0]+b*y[1]+c*y[2];u[1]=b*y[0]+a[1]*y[1]+b*y[2]+c*y[3]; //前两列和最后两列单独拿出来求中D间的循环求for(i=2;i<N-2;i++){u[i]=c*y[i-2]+b*y[i-1]+a[i]*y[i]+b*y[i+1]+c*y[i+2];}u[N-2]=c*y[N-4]+b*y[N-3]+a[N-2]*y[N-2]+b*y[N-1];u[N-1]=c*y[N-3]+b*y[N-2]+a[N-1]*y[N-1];/***求beta***/double Beta=0;for(i=0;i<N;i++){Beta+=y[i]*u[i];}//cout<<"Beta"<<k<<"="<<Beta<<"\t"; 输出每次迭代的beta /***求误差***/error=abs(Beta-Beta_)/abs(Beta);if(error<=E) //若迭代误差在精度水平内则可以停止迭代{return Beta;} //控制显示位数Beta_=Beta; //第个eta的值都要保存下来,为了与后个值进行误差计算 }if(k==2000){cout<<"error"<<endl;return 0;} //若在最大迭代次数范围内都不能满足精度要求说明不收敛}/**************反幂法求最小绝对特¬征值**************/double fmifa(){int i,k,t;double u[N],y[N]={0},yy[N]={0},b=0.16,c=-0.064,Beta_=0,error;for(i=0;i<501;i++)u[i]=1; //令u[N]=1for(k=1;k<2000;k++){double sum_u=0,gh_sum_u;for(i=0;i<N;i++){sum_u+=u[i]*u[i]; }gh_sum_u=sqrt(sum_u);for(i=0;i<N;i++){y[i]=u[i]/gh_sum_u;yy[i]=y[i]; //用重新赋值,避免求解方程组的时候改变y的值}/****LU分解法解方程组Au=y,求新的***/LU();for(i=2;i<=N;i++){double temp_b=0;for(t=max(1,i-r);t<=i-1;t++)temp_b+=cc[i-t+s][t-1]*yy[t-1];yy[i-1]=yy[i-1]-temp_b;}u[N-1]=yy[N-1]/cc[s][N-1];for(i=N-1;i>=1;i--){double temp_u=0;for(t=i+1;t<=min(i+s,N);t++)temp_u+=cc[i-t+s][t-1]*u[t-1];u[i-1]=(yy[i-1]-temp_u)/cc[s][i-1];}double Beta=0;for(i=0;i<N;i++){Beta+=y[i]*u[i];}error=abs(Beta-Beta_)/abs(Beta);if(error<=E){return (1/Beta);}Beta_=Beta;}if(k==2000){cout<<"error"<<endl;return 0;} }/**************求两数最大值的子程序**************/int max(int aa,int bb){return(aa>bb?aa:bb);}/**************求两数最小值的子程序**************/int min(int aa,int bb){return(aa<bb?aa:bb);}/**************求三数最大值的子程序**************/int max_3(int aa,int bb,int cc){ int tt;if(aa>bb)tt=aa;else tt=bb;if(tt<cc) tt=cc;return(tt);}/**************LU分解**************/void LU(){int i,j,k,t;double b=0.16,c=-0.064;/**赋值压缩后矩阵cc[5][501]**/for(i=2;i<N;i++)cc[0][i]=c;for(i=1;i<N;i++)cc[1][i]=b;for(i=0;i<N;i++)cc[2][i]=a[i];for(i=0;i<N-1;i++)cc[3][i]=b;for(i=0;i<N-2;i++)cc[4][i]=c;for(k=1;k<=N;k++){for(j=k;j<=min(k+s,N);j++){double temp=0;for(t=max_3(1,k-r,j-s);t<=k-1;t++)temp+=cc[k-t+s][t-1]*cc[t-j+s][j-1];cc[k-j+s][j-1]=cc[k-j+s][j-1]-temp;}//if(k<500){for(i=k+1;i<=min(k+r,N);i++){double temp2=0;for(t=max_3(1,i-r,k-s);t<=k-1;t++)temp2+=cc[i-t+s][t-1]*cc[t-k+s][k-1];cc[i-k+s][k-1]=(cc[i-k+s][k-1]-temp2)/cc[s][k-1];}}}}三、程序结果。

数值计算方法试题库及答案解析

数值计算方法试题库及答案解析

y 2y, y(0) 1,试问为保证该公式绝对稳定,步长 h 的取值范围为(
)。
(1) 0 h 2 , (2) 0 h 2 , (3) 0 h 2 , (4) 0 h 2
三、1、(8 分)用最小二乘法求形如 y a bx2 的经验公式拟合以下数据:
2
是否为插值型求积公式?为什么?其
代数精度是多少?
七、(9 分)设线性代数方程组 AX b 中系数矩阵 A 非奇异, X 为精确解, b 0 ,若向
~
~
量 X 是 AX b 的 一 个 近 似 解 , 残 向 量 r b A X , 证 明 估 计 式 :
~
X X
r cond ( A)
五、(8 分)已知求 a (a 0) 的迭代公式为:
1
a
xk1 2 (xk xk )
x0 0 k 0,1,2
证明:对一切 k 1,2,, xk a ,且序列xk 是单调递减的,
从而迭代过程收敛。
3 f (x)dx 3 [ f (1) f (2)]
六、(9 分)数值求积公式 0
六、(下列 2 题任选一题,4 分) 1、 1、 数值积分公式形如
1
0 xf (x)dx S(x) Af (0) Bf (1) Cf (0) Df (1)
(1) (1) 试确定参数 A, B,C, D 使公式代数精度尽量高;(2)设
1
f (x) C 4[0,1] ,推导余项公式 R(x) 0 xf (x)dx S(x) ,并估计误差。
i 1
的高斯(Gauss)型求积公式具有最高代数精确度的次
数为 2n 1。 (

数值计算大作业

数值计算大作业

数值计算大作业题目一、非线性方程求根1.题目假设人口随时间和当时人口数目成比例连续增长,在此假设下人口在短期内的增长建立数学模型。

(1)如果令()N t 表示在t 时刻的人口数目,β表示固定的人口出生率,则人口数目满足微分方程()()dN t N t dt β=,此方程的解为0()=tN t N e β; (2)如果允许移民移入且速率为恒定的v ,则微分方程变成()()dN t N t vdt β=+, 此方程的解为0()=+(1)t t vN t N e e βββ-;假设某地区初始有1000000人,在第一年有435000人移入,又假设在第一年年底该地区人口数量1564000人,试通过下面的方程确定人口出生率β,精确到410-;且通过这个数值来预测第二年年末的人口数,假设移民速度v 保持不变。

4350001564000=1000000(1)e e βββ+-2.数学原理采用牛顿迭代法,牛顿迭代法的数学原理是,对于方程0)(=x f ,如果)(x f 是线性函数,则它的求根是很容易的,牛顿迭代法实质上是一种线性化方法,其基本思想是将非线性方程0)(=x f 逐步归结为某种线性方程来求解。

设已知方程0)(=x f 有近似根k x (假定0)(≠'x f ),将函数)(x f 在点k x进行泰勒展开,有.))(()()(⋅⋅⋅+-'+≈k k k x x x f x f x f于是方程0)(=x f 可近似地表示为))(()(=-'+k k x x x f x f这是个线性方程,记其根为1k x +,则1k x +的计算公式为)()(1k k k k x f x f x x '-==+,,,2,1,0⋅⋅⋅=k这就是牛顿迭代法,简称牛顿法。

3.程序设计作出函数的图像,大概估计出根的位置fplot('1000*exp(x)+(435*x)*(exp(x)-1)-1564',[0 3]);grid大概估计出初始值x=0.5function [p1,err,k,y]=newton(f,df,p0,delta,max1) % f 是非线性系数 % df 是f 的微商 % p0是初始值% dalta 是给定允许误差 % max1是迭代的最大次数 % p1是牛顿法求得的方程近似解 % err 是p0误差估计 % k 是迭代次数 p0,feval('f',p0) for k=1:max1p1=p0-feval('f',p0)/feval('df',p0); err=abs(p1-p0); p0=p1;p1,err,k,y=feval('f',p1) if(err<delta)|(y==0), break,endp1,err,k,y=feval('f',p1) endfunction y=f(x)y=1000000*exp(x)+435000*(exp(x)-1)/x-1564000; function y=df(x)y=1000000*exp(x)+435000*(exp(x)/x-(exp(x)-1)/x^2);4.结果分析与讨论newton('f','df',1.2,10^(-4),10) 运行后得出结果 p0 =0.5000p1 =0.1679 err =0.3321 k =1 y =9.2415e+004 p1 =0.1031 err =0.0648 k =2 y =2.7701e+003 p1 =0.1010 err =0.0021 k =3 y =2.6953p1 =0.1010 err =2.0129e-006 k =4 y = 2.5576e-006 ans =0.1010运算后的结果为1010.0=β,通过这个数值来预测第二年年末的人口数,0.10100.1010435000f(t)=1000000(1)0.1010t te e +-t=2时候对于f ()2187945.865x =实践表明,当初始值难以确定时,迭代法就不一定收敛了,因此要根据问题实际背景或者二分法先得一个较好的初始值,然后再进行迭代;再者迭代函数选择不合适的话,采用不动点迭代法也有可能出现不收敛的情况;因此我采用的是牛顿法。

数值计算方法大作业

数值计算方法大作业

题目利用数值计算方法求取基尼系数姓名与学号指导教师年级与专业所在学院一、问题综述:基尼系数(Gini coefficient),是20世纪初意大利学者科拉多·吉尼根据劳伦茨曲线所定义的判断收入分配公平程度的指标。

是比例数值,在0和1之间。

基尼指数(Gini index)是指基尼系数乘100倍作百分比表示。

在民众收入中,如基尼系数最大为“1”,最小等于“0”。

前者表示居民之间的收入分配绝对不平均(即所有收入都集中在一个人手里,其余的国民没有收入),而后者则表示居民之间的收入分配绝对平均,即人与人之间收入绝对平等,但这两种情况只出现在理论上;因此,基尼系数的实际数值只能介于0~1之间,基尼系数越小收入分配越平均,基尼系数越大收入分配越不平均。

设右图中的实际收入分配曲线(红线)和收入分配绝对平等线(绿线)之间的面积为A,和收入分配绝对不平等线(蓝线)之间的面积为B,则表示收入与人口之间的比例的基尼系数为AA+B。

如果A为零,即基尼系数为0,表示收入分配完全平等(红线和绿线重叠);如果B为零,则系数为1,收入分配绝对不平等(红线和蓝线重叠)。

该系数可在0和1之间取任何值。

实际上,一般国家的收入分配,既不是完全平等,也不是完全不平等,而是在两者之间,劳伦茨曲线为一条凸向横轴的曲线。

收入分配越趋向平等,劳伦茨曲线的弧度越小(斜度越倾向45度),基尼系数也越小;反之,收入分配越趋向不平等,劳伦茨曲线的弧度越大,那么基尼系数也越大。

基尼系数的调节需要国家通过财政政策进行国民收入的二次分配,例如对民众的财政公共服务支出和税收等,从而让收入均等化,令基尼系数缩小。

基尼系数由于给出了反映居民之间贫富差异程度的数量界线,可以较客观、直观地反映和监测居民之间的贫富差距,预报、预警和防止居民之间出现贫富两极分化。

因此得到世界各国的广泛认同和普遍采用。

联合国有关组织规定:●若低于0.2表示收入平均;●0.2-0.3表示相对平均;●0.3-0.4表示相对合理;●0.4-0.5表示收入差距大;●0.6以上表示收入差距悬殊。

数值计算方法答案

数值计算方法答案

数值计算方法习题一(2)习题二(6)习题三(15)习题四(29)习题五(37)习题六(62)习题七(70)2009.9,9习题一1.设x >0相对误差为2%,4x 的相对误差。

解:由自变量的误差对函数值引起误差的公式:(())(())'()()()()f x xf x f x x f x f x δδ∆=≈得(1)()f x =11()()*2%1%22x x δδδ≈===;(2)4()f x x =时444()()'()4()4*2%8%x x x x x x δδδ≈===2.设下面各数都是经过四舍五入得到的近似数,即误差不超过最后一位的半个单位,试指出他们各有几位有效数字。

(1)12.1x =;(2)12.10x =;(3)12.100x =。

解:由教材9P 关于1212.m nx a a a bb b =±型数的有效数字的结论,易得上面三个数的有效数字位数分别为:3,4,53.用十进制四位浮点数计算 (1)31.97+2.456+0.1352; (2)31.97+(2.456+0.1352)哪个较精确?解:(1)31.97+2.456+0.1352 ≈21((0.3197100.245610)0.1352)fl fl ⨯+⨯+ =2(0.3443100.1352)fl ⨯+=0.3457210⨯(2)31.97+(2.456+0.1352)21(0.319710(0.245610))fl fl ≈⨯+⨯ = 21(0.3197100.259110)fl ⨯+⨯ =0.3456210⨯易见31.97+2.456+0.1352=0.345612210⨯,故(2)的计算结果较精确。

4.计算正方形面积时,若要求面积的允许相对误差为1%,测量边长所允许的相对误差限为多少?解:设该正方形的边长为x ,面积为2()f x x =,由(())(())'()()()()f x xf x f x x f x f x δδ∆=≈解得(())()()'()f x f x x xf x δδ≈=2(())(())22f x x f x x xδδ==0.5%5.下面计算y 的公式哪个算得准确些?为什么?(1)已知1x <<,(A )11121xy x x-=-++,(B )22(12)(1)x y x x =++; (2)已知1x >>,(A)y =,(B)y = (3)已知1x <<,(A )22sin x y x =,(B )1cos 2xy x-=;(4)(A)9y =(B)y =解:当两个同(异)号相近数相减(加)时,相对误差可能很大,会严重丧失有效数字;当两个数相乘(除)时,大因子(小除数)可能使积(商)的绝对值误差增大许多。

数值计算方法大作业--资料

数值计算方法大作业--资料

数值计算方法大作业--资料-CAL-FENGHAI-(2020YEAR-YICAI)_JINGBIAN计算方法大作业学生学号: ********学生姓名: ****专业班级: ***********摘要:大作业通过MATLAB在计算方法中的应用实例,探讨了MATLAB在计算方法中的应用方法和技巧,对运用计算机软件完成“计算方法”课程的图形绘制,多项式方程的求解,计算方法分析具有较好的参考价值。

关键字:MATLAB应用迭代法多项式引言在科学研究与工程设计中,经常会遇到数学模型的求解问题,然而在许多情况下,要获得模型问题的准确解是十分困难的,甚至是不可能的。

因此,研究各种数学问题的近似解法非常重要。

数值计算方法又称计算方法或数值计算分析,是一门与计算机应用密切结合的实用性很强的数学课程。

数值计算方法提供的算法具有以下特点:1.面向计算机,根据计算机的特点设计可行的算法。

2.有可靠的理论依据。

3.高效率。

数值计算方法既重视与方法有关的理论,又重视方法的实际运用,而且数值计算方法课程涉及的面较广泛,包括了微积分、线性代数、常微积分方程等数学问题的数值方法。

所以我们只有努力的掌握这几门课程的基本内容,才能学好这门课程。

掌握数值计算方法,包括数组和数组函数,矩阵和矩阵函数的创建与操作,关系与逻辑操作符的运算,多项式计算,数据分析,以及方程与方程组的解法。

掌握Matla图形和3D可视化的技术,围绕数据成图机理,绘图要旨和修饰技法熟悉各种绘图指令和交互操作工具。

包括二维,三维和高维图形绘制,图形的色彩,光源和材质等效果的处理,以及图形句柄操作和动画制作技术。

Matlab数值计算,数值计算功能是Matlab最具代表性的特点,也是最基本、最重要的功能,它是备受欢迎的基石。

Matlab能够成为世界上最优秀的数学软件之一和它出色的数值运算能力是分不开的。

Matlab在数值运算中以数组和矩阵为基础。

数组是Matlab运算中一个重要的数据组织形式。

数值计算方法习题(2018-解答仅供参考)

数值计算方法习题(2018-解答仅供参考)

习 题 一1、下列各数都是经过四舍五入得到的近似值。

试分别指出它们的绝对误差、相对误差和有效数字的位数。

35801=x ,00476.02=x ,33101430.0⨯=x ,24102958-⨯=x ,85000.55=x 。

解:2、已知2031.1=a ,978.0=b 是经过四舍五入后得到的近似值,问b a b a ⨯+,有几位有效数字?解:321110,1022a ab b *-*--≤⨯-≤⨯,而 2.1811, 1.1766a b a b +=⨯=()()3212111101010222a b a b a a b b ****---+-+≤-+-≤⨯+⨯≤⨯故a b +至少具有2位有效数字。

()()32120.978 1.2031110100.006510222ab a b b a a a b b *****----≤-+-≤⨯+⨯=≤⨯ 故a b +至少具有2位有效数字。

3、求二次方程01162=+-x x 的较小正根,取94.763≈,要求有3位有效数字。

解:*1228887.940.06x x x ==-==,*2x 只有一位有效数字。

若改用2180.062715.94x =≈,具有三位有效数字。

4、正方形的边长约cm 100,问测量边长时误差应多大,才能保证面积的误差不超过21cm ?解:正方形的面积函数为2()A x x =(*)2*(*)A A x εε∴=.当*100x =时,若(*)1A ε≤, 则21(*)102x ε-≤⨯故测量中边长误差限不超过0.005cm 时,才能使其面积误差不超过21cm 5、改变下列表达式,使其计算结果比较精确。

(1)1,11>>--+x xx x x ; (2)1),1ln(2>>--x x x ; (3)1,1<<-x e x ; (4)xxsin cos 1-。

解:(1(2) (3)(4)22sin 1cos 2tan sin 22sin cos 22x x xx x x -==习 题 二1、已知314567.032.0sin =,333487.034.0sin =,352274.036.0sin =,用线性插值及抛物插值计算3367.0sin 的值并估计截断误差。

数值计算方法丁丽娟课后习题答案

数值计算方法丁丽娟课后习题答案

数值计算方法丁丽娟课后习题答案【篇一:北京理工大学数值计算方法大作业数值实验1】)书p14/4分别将区间[?10,10]分为100,200,400等份,利用mesh或surf命令画出二元函数的三维图形。

z=|??|+ ??+?? +??++??【matlab求解】[x,y]=meshgrid(-10:0.1:10);a=exp(-abs(x));b=cos(x+y);c=1./(x.^2+y.^2+1);z=a+b+c;mesh(x,y,z);[x,y]=meshgrid(-10:0.05:10);a=exp(-abs(x));b=cos(x+y);c=1./(x.^2+y.^2+1);z=a+b+c;mesh(x,y,z);[x,y]=meshgrid(-10:0.025:10); a=exp(-abs(x));b=cos(x+y);c=1./(x.^2+y.^2+1);z=a+b+c;mesh(x,y,z);(二)书p7/1.3.2数值计算的稳定性(i)取= ??c语言程序—不稳定解 +=ln1.2,按公式=?? (n=1,2,…) #includestdio.h#includeconio.h#includemath.hvoid main(){float m=log(6.0)-log(5.0),n;int i;i=1;printf(y[0]=%-20f,m); while(i20){n=1/i-5*m;printf(y[%d]=%-20f,i,n);m=n;i++;if (i%3==0) printf(\n); }getch();}(ii) c语言程序—稳定解≈??[ ??+?? +?? ??+??按公式 =??(??)#includestdio.h#includeconio.h#includemath.hvoid main(){float m=(1/105.0+1/126.0)/2,n; k=n,n-1,n-2,…)(【篇二:北京理工大学数值计算方法大作业数值实验4】 p260/1考纽螺线的形状像钟表的发条,也称回旋曲线,它在直角坐标系中的参数方程为= ?????????????????? ?? ??????????= ?????????????? ??曲线关于原点对称,取a=1,参数s的变化范围[-5,5],容许误差限分别是,,和。

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

目录第一章非线性方程求根 (3)1.1迭代法 (3)1.2牛顿法 (4)1.3弦截法 (5)1.4二分法 (6)第二章插值 (7)2.1线性插值 (7)2.2二次插值 (8)2.3拉格朗日插值 (9)2.4分段线性插值 (10)2.5分段二次插值 (11)第三章数值积分 (13)3.1复化矩形积分法 (13)3.2复化梯形积分法 (14)3.3辛普森积分法 (15)3.4变步长梯形积分法 (16)第四章线性方程组数值法 (17)4.1约当消去法 (17)4.2高斯消去法 (18)4.3三角分解法 (20)4.4雅可比迭代法 (21)4.5高斯—赛德尔迭代法 (23)第五章常积分方程数值法 (25)5.1显示欧拉公式法 (25)5.2欧拉公式预测校正法 (26)5.3改进欧拉公式法 (27)5.4四阶龙格—库塔法 (28)数值计算方法第一章非线性方程求根1.1迭代法程序代码:Private Sub Command1_Click()x0 = Val(InputBox("请输入初始值x0"))ep = Val(InputBox(请输入误差限ep))f = 0While f = 0X1 = (Exp(2 * x0) - x0) / 5If Abs(X1 - x0) < ep ThenPrint X1f = 1Elsex0 = X1End IfWendEnd Sub例:求f(x)=e2x-6x=0在x=0.5附近的根(ep=10-10)1.2牛顿法程序代码:Private Sub Command1_Click()b = Val(InputBox("请输入被开方数x0"))ep = Val(InputBox(请输入误差限ep))f = 0While f = 0X1 = x0 - (x0 ^ 2 - b) / (2 * b)If Abs(X1 - x0) < ep ThenPrint X1f = 1Elsex0 = X1End IfWendEnd Sub例:求56的值。

(ep=10-10)1.3弦截法程序代码:Private Sub Command1_Click()x0 = Val(InputBox("请输入第一个初始值x0"))X1 = Val(InputBox("请输入第二个初始值x1"))ep = Val(InputBox("请输入误差限ep"))f = 0While f = 0X2 = X1 - (X1 ^ 8 - 13) * (X1 - x0) / ((X1 ^ 8 - 13) - (x0 ^ 8 - 13))If Abs(X2 - X1) < ep ThenPrint X2f = 1Elsex0 = X1X1 = X2End IfWendEnd Sub例:求f(x)=x8-13的正根(初始值x1=1,x2=10,ep=10-10)1.4二分法程序代码:Private Sub Command1_Click()a = Val(InputBox("请输入区间端点a"))b = Val(InputBox("请输入区间端点b"))ep = Val(InputBox("请输入误差限ep"))f = 0While f = 0x = (a + b) / 2fx = Exp(-x / 7) * (9 - 2 * x) - 8fa = Exp(-a / 7) * (9 - 2 * a) - 8If fx = 0 Thenf = 1Print "方程的根是", xElseIf fa * fx > 0 Thena = xElseb = xEnd IfIf Abs(b - a) < ep Thenx = (b + a) / 2f = 1Print "方程的根是", xEnd IfEnd IfWendEnd Sub例:求方程f(x)=e-7/x(9-2x)-8在区间[0,1]的实根。

(ep=10-10)第二章插值2.1线性插值程序代码:Private Sub Command1_Click()X0 = Val(InputBox("请输入第一个结点X:"))Y0 = Val(InputBox("请输入第一个结点Y:"))X1 = Val(InputBox("请输入第二个结点X:"))Y1 = Val(InputBox("请输入第二个结点Y:"))f = 0While f = 0x = Val(InputBox("请输入未知点的自变量值X:"))L0 = (x - X1) / (X0 - X1)L1 = (x - X0) / (X1 - X0)y = L0 * Y0 + L1 * Y1Print "x="; x, "y="; yf = Val(InputBox("是否继续(0/1):"))WendEnd Sub例:已知两点(13 , 1)、(49 , 8),求30处的值。

2.2二次插值程序代码:Private Sub Command1_Click()X0 = Val(InputBox("请输入第一个结点X:"))Y0 = Val(InputBox("请输入第一个结点Y:"))X1 = Val(InputBox("请输入第二个结点X:"))Y1 = Val(InputBox("请输入第二个结点Y:"))X2 = Val(InputBox("请输入第三个结点X:"))Y2 = Val(InputBox("请输入第三个结点Y:"))f = 0While f = 0x = Val(InputBox("请输入未知点的自变量值X:"))L0 = (x - X1) * (x - X2) / (X0 - X1) / (X0 - X2)L1 = (x - X0) * (x - X2) / (X1 - X0) / (X1 - X2)L2 = (x - X0) * (x - X1) / (X2 - X0) / (X2 - X1)y = L0 * Y0 + L1 * Y1 + L2 * Y2Print "x="; x, "y="; yf = Val(InputBox("是否继续(0/1):"))WendEnd Sub例:已知三点(81 ,9)、(100 ,10)、(121 ,10),求98处的值。

2.3拉格朗日插值程序代码:Private Sub Command1_Click()Dim x(), y()n = Val(InputBox("请输入插值节点数N"))ReDim x(n), y(n)For i = 0 To nx(i) = Val(InputBox("请输入插值节点x(" + Str(i) + ")"))y(i) = Val(InputBox("请输入插值节点y(" + Str(i) + ")"))Next if = 0While f = 0xx = Val(InputBox("请输入未知点的自变量x:"))Sum = 0For i = 0 To nt = 1For j = 0 To nIf j <> i Thent = t * (xx - x(j)) / (x(i) - x(j))End IfNext jSum = Sum + t * y(i)Next iPrint "x="; xx, "y="; Sumf = Val(InputBox("是否继续(0/1)"))WendEnd Sub例:已知四点(100 ,10)、(81 ,9)、(64 ,8)、(49 ,7),求87 处的值。

2.4分段线性插值程序代码:Private Sub Command1_Click()Dim x(), y()n = Val(InputBox("请输入插值节点数N"))ReDim x(n), y(n)For i = 0 To nx(i) = Val(InputBox("请输入插值节点x(" + Str(i) + ")")) y(i) = Val(InputBox("请输入插值节点y(" + Str(i) + ")")) Next if = 0While f = 0xx = Val(InputBox("请输入未知点的自变量x:"))L = 0j = 1While L = 0If xx < x(j) Thenk = j + 1L = 1Elsej = j + 1If j > n - 1 Thenk = n - 1L = 1End IfEnd IfWendl0 = (xx - x(k)) / (x(k - 1) - x(k))l1 = (xx - x(k - 1)) / (x(k) - x(k - 1))yy = l0 * y(k - 1) + l1 * y(k)Print "x="; xx, "y="; yyf = Val(InputBox("是否继续(0/1)"))WendEnd Sub例:已知三点(361 , 19)、(324 ,18)、(289 ,17),N=2,求300处的值。

2.5分段二次插值程序代码:Private Sub Command1_Click()Dim x(), y()n = Val(InputBox("请输入插值节点数N"))ReDim x(n), y(n)For i = 0 To nx(i) = Val(InputBox("请输入插值节点x(" + Str(i) + ")"))y(i) = Val(InputBox("请输入插值节点y(" + Str(i) + ")"))Next if = 0While f = 0xx = Val(InputBox("请输入未知点的自变量x:"))If x0 < x(1) Thenk = 1f = 1End Ifi = 2Do While f = 0 And i >= n - 1If x0 < x(i) ThenIf x0 - x(i - 1) < x(i) - x0 Thenk = i - 1f = 1Elsek = if = 1End IfElsei = i + 1End IfLoopIf f = 0 Thenk = n - 1End Ifl1 = (xx - x(k + 1)) * (xx - x(k)) / ((x(k - 1) - x(k + 1)) * (x(k - 1) - x(k)))l2 = (xx - x(k + 1)) * (xx - x(k - 1)) / ((x(k) - x(k + 1)) * (x(k) - x(k - 1)))l3 = (xx - x(k)) * (xx - x(k - 1)) / ((x(k + 1) - x(k)) * (x(k + 1) - x(k - 1)))yy = l1 * y(k - 1) + l2 * y(k) + l3 * y(k + 1)Print "x="; xx, "y="; yyf = Val(InputBox("是否继续(0/1)"))WendEnd Sub例:已知三点(225 , 15)、(196 ,14)、(169 ,13),求180处的值。

相关文档
最新文档