数值计算方法上机实习题
数值计算方法上机实习题NEW
数值计算方法上机实习题1. 设⎰+=105dx x x I nn , (1) 由递推公式nI I n n 151+-=-,从0I 的几个近似值出发,计算20I ; (2) 粗糙估计20I ,用n I I n n 51511+-=-,计算0I ; (3) 分析结果的可靠性及产生此现象的原因(重点分析原因)。
2. 求方程0210=-+x e x 的近似根,要求41105-+⨯<-k k x x ,并比较计算量。
(1) 在[0,1]上用二分法;(2) 取初值00=x ,并用迭代1021x k e x -=+; (3) 加速迭代的结果;(4) 取初值00=x ,并用牛顿迭代法;(5) 分析绝对误差。
3.钢水包使用次数多以后,钢包的容积增大,数据如下:(注:增速减少,用何种模型) 4.设⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛----------------=410100141010014101101410010141001014A ,⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛--=625250b ,b x =A 分析下列迭代法的收敛性,并求42110-+≤-k k x x 的近似解及相应的迭代次数。
(1) JACOBI 迭代;(2) GAUSS-SEIDEL 迭代;(3) SOR 迭代(95.0,95.1,334.1=ω)。
5.用逆幂迭代法求⎪⎪⎪⎭⎫ ⎝⎛=111123136A 最接近于11的特征值和特征向量,准确到310-。
6.用经典R-K 方法求解初值问题(1)⎩⎨⎧-+-='++-='x x y y y x y y y sin 2cos 22sin 22212211,]10,0[∈x , ⎩⎨⎧==3)0(2)0(21y y ; (2)⎩⎨⎧-+-='++-='x x y y y x y y y sin 999cos 999999998sin 22212211,]10,0[∈x , ⎩⎨⎧==3)0(2)0(21y y 。
数值计算方法上机实习题答案
1. 设⎰+=105dx x x I nn , (1) 由递推公式n I I n n 151+-=-,从0I 的几个近似值出发,计算20I ; 解:易得:0I =ln6-ln5=0.1823,程序为:I=0.182;for n=1:20I=(-5)*I+1/n;endI输出结果为:20I = -3.0666e+010(2) 粗糙估计20I ,用n I I n n 515111+-=--,计算0I ; 因为 0095.056 0079.010********≈<<≈⎰⎰dx x I dx x 所以取0087.0)0095.00079.0(2120=+=I 程序为:I=0.0087;for n=1:20I=(-1/5)*I+1/(5*n);endI0I = 0.0083(3) 分析结果的可靠性及产生此现象的原因(重点分析原因)。
首先分析两种递推式的误差;设第一递推式中开始时的误差为000I I E '-=,递推过程的舍入误差不计。
并记nn n I I E '-=,则有01)5(5E E E n n n -==-=- 。
因为=20E 20020)5(I E >>-,所此递推式不可靠。
而在第二种递推式中n n E E E )51(5110-==-= ,误差在缩小,所以此递推式是可靠的。
出现以上运行结果的主要原因是在构造递推式过程中,考虑误差是否得到控制,即算法是否数值稳定。
2. 求方程0210=-+x e x 的近似根,要求41105-+⨯<-k k x x ,并比较计算量。
(1) 在[0,1]上用二分法;程序:a=0;b=1.0;while abs(b-a)>5*1e-4if exp(c)+10*c-2>0b=c;else a=c;endendc结果:c =0.0903(2) 取初值00=x ,并用迭代1021x k e x -=+; 程序:x=0;a=1;while abs(x-a)>5*1e-4a=x;x=(2-exp(x))/10;endx结果:x =0.0905(3) 加速迭代的结果;程序:x=0;a=0;b=1;while abs(b-a)>5*1e-4a=x;y=exp(x)+10*x-2;z=exp(y)+10*y-2;x=x-(y-x)^2/(z-2*y+x);b=x;endx结果:x =0.0995(4) 取初值00=x ,并用牛顿迭代法;程序:x=0;a=0;b=1;while abs(b-a)>5*1e-4x=x-(exp(x)+10*x-2)/(exp(x)+10);b=x;endx结果:x =0.0905(5) 分析绝对误差。
数值计算上机实习题目(matlab编程)
数值计算上机实习题目(matlab编程)非线性方程求根一、实验目的本次实验通过上机实习,了解迭代法求解非线性方程数值解的过程和步骤。
二、实验要求1、用迭代法求方程230x x e -=的根。
要求:确定迭代函数?(x),使得x=?(x),并求一根。
提示:构造迭代函数2ln(3)x ?=。
2、对上面的方程用牛顿迭代计算。
3、用割线法求方程3()310f x x x =--=在02x =附近的根。
误差限为410-,取012, 1.9x x ==。
三、实验内容1、(1)首先编写迭代函数,记为iterate.mfunction y=iterate(x)x1=g(x); % x 为初始值。
n=1;while(abs(x1-x)>=1.0e-6)&(n<=1000) % 迭代终止的原则。
x=x1;x1=g(x);n=n+1;endx1 %近似根n %迭代步数(2)后编制函数文件?(x),记为g.mfunction y=g(x)y=log(3*x.^2);(3)设初始值为0、3、-3、1000,观察初始值对求解的影响。
将结果记录在文档中。
>>iterate(0)>>iterate(3) 等等2、(1)首先编制牛顿迭代函数如下,记为newton.mfunction y=newton(x0)x1=x0-fc(x0)/df(x0); % 牛顿迭代格式n=1;while(abs(x1-x0)>=1.0e-6)&(n<=1000000) % 迭代终止的原则。
x0=x1;x1=x0-fc(x0)/df(x0);n=n+1;endx1 %近似根n %迭代步数(2)对题目中的方程编制函数文件,记为fc.mfunction y=fc(x)y=3*x.^2-exp(x)编制函数的导数文件,记为df.mfunction y=df(x)y=6*x-exp(x)(3)在MATLAB 命令窗计算,当设初始值为0时,newton(0);给定不同的初始值,观察用牛顿法求解时所需要的迭代步数,并与上面第一题的迭代步数比较。
数值计算方法上机题目资料
2. 程序输入、输出用文件形式。 3. 编程语言要求用C, 编程环境TC或VC 4. 程序要求调试通过。 5. 每个方法要求给出一个具体的算例(可选
作业题)来验证。
五、上机报告要求
1.报告内容包括:
每种方法的算法原理及程序框图。 程序使用说明 具体算例及结果
上机调试体会及收获。
2.报告要手写。
六、上机报告及源程序提交时间
1.上机报告在考试当天提交。 2.源程序在考试前提交。
提交格式:文件夹(班级+姓名)
输入文件 程序文件夹 输出文件
不要拷贝其它文件!!! 源程序
六、上机报告及源程序提交时间
源程序提交: 把以上文件压缩后,发送到以下邮箱:
haoyq@
七、考核方式
1.算法手算笔试(80%)+上机内容笔试 (10%)+上机报告(10%)
2.上机内容笔试可能形式:
编一段算法程序 给出一段算法程序,说明算法的名称。 程序填空 程序改错(包括算法和语法的错误)
数值计算方法上机练习
一、上机练习目的
复习和巩固数值计算方法的基本数学模型, 全面掌握运用计算机进行数值计算的具体 过程及相关问题。
利用计算机语言独立编写、调试数值计算 方法程序,培养学生利用计算机和所学理 论知识分析解决实际问题的能力。
二、上机练习任务
• 利用计算机基本C语言编写并调试一系列 数值方法计算通用程序,并能正确计算给 定题目,掌握调试技能。
• 掌握文件使用编程技能,如文件的各类操 作,数据格式设计、通用程序运行过程中 文件输入输出运行方式设计等。
• 写出上机练习报告。
三、数值计算方法上机题目
数值计算方法上机实验考试答案
数值计算方法上机实验考试答案1. 小型火箭初始质量为900千克,其中包括600千克燃料。
火箭竖直向上发射时燃料以15千克/秒的速率燃烧掉,由此产生30000牛顿的恒定推力。
当燃料用尽时引擎关闭。
设火箭上升的整个过程中,空气阻力与速度平方成正比,比例系数为0.4(千克/米)。
重力加速度取9.8米/秒2.A. 建立火箭升空过程的数学模型(微分方程);B. 求引擎关闭瞬间火箭的高度、速度、加速度,及火箭到达最高点的时间和高度。
解:A . 建立模型:)(t x ——t 时刻的火箭高度;T =30000(牛顿)——火箭推力,当时间t >40秒时,T=0; m t m 150-=——火箭飞行过程中的质量,t >40秒时,300=m 千克0m =900(千克)——火箭初始质量; 1m =600(千克)——燃料质量;c =15(千克/秒)——燃料的燃烧速率; k =0.4(千克/米)——空气阻力系数; g =9.8(米/秒2)——重力加速度由能量守恒定律,可得到火箭飞行过程的方程:mg T t x k t x m -+'-=''2)]([)(这是一个初值问题,初始条件为 0)0(,0)0(='=x x设)(),(21t x x t x x '==,则问题化为求下列微分方程组的初值问题:⎪⎩⎪⎨⎧--+--='='g t m T x t m k x x x 151******** 0)0(,0)0(21==x xB . 关闭引擎时4015/600==t (秒),所求的是此时火箭的高度)40(1x ;速度)40(2x ;加速度)40()40(21x x '''或,及火箭到最高点的时间m t 和高度)(1m t x 。
具体的Matlab 程序如下: 首先建立微分方程的的m 文件:function y=huojian(t,x)k=0.4;g=9.8;m0=900;T=30000;m=m0-15*t;if t>40T=0;m=300;endy=[x(2),-(k/m)*x(2)^2+T/m-g]';主程序:%feixing.mk=0.4;g=9.8;m0=900;T=30000;x0=[0,0];ts=0:1:55;[t,x]=ode23('huojian',ts,x0);[t,x(:,1)]%------a=[t,x];x40=a(41,2) %燃料用尽时的高度v40=a(41,3) %燃料用尽时的速度a40=-(k/300)*v40^2+T/300-g %燃料用尽时的加速度%-------xmax=max(x(:,1)) %火箭到达最高点的高度subplot(2,1,1),plot(t,x(:,1)),title('altitude') subplot(2,1,2),plot(t,x(:,2)),title('speed')运行结果为:1.0e+003 *0 00.0010 0.01180.0020 0.04750.0030 0.10670.0040 0.18890.0050 0.29270.0060 0.41680.0070 0.55920.0080 0.71800.0090 0.89140.0100 1.07700.0380 7.80510.0390 8.06310.0400 8.32240.0410 8.54000.0420 8.70040.0430 8.82420.0440 8.92180.0450 8.99940.0460 9.06070.0470 9.10830.0480 9.14370.0490 9.16810.0500 9.18220.0510 9.18640.0520 9.18070.0530 9.16510.0540 9.13900.0550 9.1018x40 =8322.4v40 = 254.1728a40 = 4.0616xmax =9186.4关闭引擎时4015/600==t (秒),此时火箭的高度h=)40(1x =8322.4米,速度v=)40(2x =254.1728米/秒,加速度为a=)40(1x ''= 4.0616米/秒2,火箭到最高点的时间m t =51秒,高度)(1m t x =9186.4米。
数值分析计算实习题答案
数值分析计算实习题答案数值分析计算实习题答案数值分析是一门研究如何利用计算机对数学问题进行近似求解的学科。
在数值分析的学习过程中,实习题是一种重要的学习方式,通过实践来巩固理论知识,并培养解决实际问题的能力。
本文将为大家提供一些数值分析计算实习题的答案,希望能够帮助大家更好地理解和掌握数值分析的相关知识。
一、插值与拟合1. 已知一组数据点,要求通过这些数据点构造一个一次插值多项式,并求出在某一特定点的函数值。
答案:首先,我们可以根据给定的数据点构造一个一次插值多项式。
假设给定的数据点为(x0, y0), (x1, y1),我们可以构造一个一次多项式p(x) = a0 + a1x,其中a0和a1为待定系数。
根据插值条件,我们有p(x0) = y0,p(x1) = y1。
将这两个条件代入多项式中,可以得到一个方程组,通过求解这个方程组,我们就可以确定a0和a1的值。
最后,将求得的多项式代入到某一特定点,就可以得到该点的函数值。
2. 已知一组数据点,要求通过这些数据点进行最小二乘拟合,并求出拟合曲线的表达式。
答案:最小二乘拟合是一种通过最小化误差平方和来找到最佳拟合曲线的方法。
假设给定的数据点为(x0, y0), (x1, y1),我们可以构造一个拟合曲线的表达式y =a0 + a1x + a2x^2 + ... + anx^n,其中a0, a1, ..., an为待定系数。
根据最小二乘拟合原理,我们需要最小化误差平方和E = Σ(yi - f(xi))^2,其中yi为实际数据点的y值,f(xi)为拟合曲线在xi处的函数值。
通过求解这个最小化问题,我们就可以确定拟合曲线的表达式。
二、数值积分1. 已知一个函数的表达式,要求通过数值积分的方法计算函数在某一区间上的定积分值。
答案:数值积分是一种通过将定积分转化为数值求和来近似计算的方法。
假设给定的函数表达式为f(x),我们可以将定积分∫[a, b]f(x)dx近似为Σwi * f(xi),其中wi为权重系数,xi为待定节点。
数值分析上机实验题参考
数值分析论文数值积分 一、问题提出选用复合梯形公式,复合Simpson 公式,Romberg 算法,计算I = dx x ⎰-4102sin 4 ()5343916.1≈II =dx x x ⎰1sin ()9460831.0,1)0(≈=I fI = dx xe x⎰+1024 I =()dx x x ⎰++1211ln 二、要求编制数值积分算法的程序;分别用两种算法计算同一个积分,并比较其结果;分别取不同步长()/ a b h -=n ,试比较计算结果(如n = 10, 20等); ﹡给定精度要求ε,试用变步长算法,确定最佳步长﹡。
三、目的和意义深刻认识数值积分法的意义; 明确数值积分精度与步长的关系;根据定积分的计算方法,可以考虑二重积分的计算问题引言一、数值求积的基本思想实际问题当中常常需要计算积分,有些数值方法。
如微分方程和积分方程的求解,也都和积分计算相联系。
依据人们熟悉的微积分基本原理,对于积分I=⎰a b f(x)dx,只要找到被积函数f(x)和原函数F(x),便有下列牛顿-莱布尼茨公式:I=⎰a b f(x)dx=F(b)-F(a).但实际使用这种求积方法往往有困难,因为大量的被积函数,诸如x xsin,2xe-等,其原函数不能用初等函数表达,故不能用上述公式计算。
另外,当f(x)是由测量或数值计算给出的一张数据表时,牛顿-莱布尼茨公式也不能直接运用,因此有必要研究积分的数值计算问题。
二、数值积分代数精度数值求积方法是近似方法,为要保证精度,我们自然希望求积公式能对“尽可能多”的函数准确成立,就提出了所谓代数精度的概念。
如果某个求积公式对次数不超过m的多项式均能准确成立,但对m+1次多项式就不能准确成立,则称该求积公式具有m次代数精度。
三、复合求积公式为了提高精度,通常可以把积分区间分成若干子区间(通常是等分),再在每个子区间用低阶求积公式,即复化求积法,比如复化梯形公式与复化辛普森公式。
计算方法上机实习题大作业(实验报告)
计算方法实验报告班级: 学号: 姓名: 成绩:1 舍入误差及稳定性一、实验目的(1)通过上机编程,复习巩固以前所学程序设计语言及上机操作指令;(2)通过上机计算,了解舍入误差所引起的数值不稳定性二、实验内容1、用两种不同的顺序计算1000021n n -=∑,分析其误差的变化 2、已知连分数()101223//(.../)n n a f b b a b a a b =++++,利用下面的算法计算f : 11,i n n i i i a d b d b d ++==+ (1,2,...,0)i n n =-- 0f d = 写一程序,读入011,,,...,,,...,,n n n b b b a a 计算并打印f3、给出一个有效的算法和一个无效的算法计算积分1041nn x y dx x =+⎰ (0,1,...,10)n = 4、设2211N N j S j ==-∑,已知其精确值为1311221N N ⎛⎫-- ⎪+⎝⎭(1)编制按从大到小的顺序计算N S 的程序(2)编制按从小到大的顺序计算N S 的程序(3)按两种顺序分别计算10001000030000,,,S S S 并指出有效位数三、实验步骤、程序设计、实验结果及分析1、用两种不同的顺序计算1000021n n -=∑,分析其误差的变化 (1)实验步骤:分别从1~10000和从10000~1两种顺序进行计算,应包含的头文件有stdio.h 和math.h(2)程序设计:a.顺序计算#include<stdio.h>#include<math.h>void main(){double sum=0;int n=1;while(1){sum=sum+(1/pow(n,2)); if(n%1000==0)printf("sun[%d]=%-30f",n,sum);if(n>=10000)break;n++;}printf("sum[%d]=%f\n",n,sum); }b.逆序计算#include<stdio.h>#include<math.h>void main(){double sum=0;int n=10000;while(1){sum=sum+(1/pow(n,2));if(n%1000==0)printf("sum[%d]=%-30f",n,sum);if(n<=1)break;n--;}printf("sum[%d]=%f\n",n,sum);}(3)实验结果及分析:程序运行结果:a.顺序计算b.逆序计算结果分析:两种不同顺序计算结果是一样的,顺序计算误差从一开始就很小,而逆序计算误差最开始十分大,后来结果正确。
西南交通大学数值分析上机实习
目录解题: (1)题目一: (1)1.1计算结果 (1)1.2结果分析 (1)题目二: (2)2.1计算结果 (2)2.2结果分析 (3)题目三: (4)3.1计算结果 (4)3.2结果分析 (5)总结 (5)附录 (6)Matlab程序: (6)题目一: (6)第一问Newton法: (6)第二问Newton法: (6)第一问Steffensen加速法: (7)第二问Steffensen加速法: (7)题目二 (8)1、Jacobi迭代法 (8)2、Causs-Seidel迭代法 (8)题目三: (9)题目一:分别用牛顿法,及基于牛顿算法下的Steffensen 加速法(1)求ln(x +sin x )=0的根。
初值x0分别取0.1, 1,1.5, 2, 4进行计算。
(2)求sin x =0的根。
初值x0分别取1,1.4,1.6, 1.8,3进行计算。
分析其中遇到的现象与问题。
1.1计算结果求ln(x +sin x )=0的根,可变行为求解x-sinx-1=0的根。
1.2结果分析从结果对比我们可发现牛顿—Steffensen 加速法比牛顿法要收敛的快,牛顿法对于初值的选取特别重要,比如第(1)问中的初值为4的情况,100次内没有迭代出来收敛解,而用Steffensen 加速法,7次迭代可得;在第(2)问中的初值为1.6的情况,收敛解得31.4159,分析其原因应该是x x f cos )('=,x0=1.62π≈,0)('≈x f ;迭代式在迭代过程中会出现分母趋近于0,程序自动停止迭代的情况,此时得到的x 往往非常大,而在第一问中我们如果转化为用x+sinx=1,则可以收敛到结果。
用雅格比法与高斯-赛德尔迭代法解下列方程组Ax=b,研究其收敛性,上机验证理论分析是否正确,比较它们的收敛速度,观察右端项对迭代收敛有无影响。
(1)A行分别为A1=[6,2,-1],A2=[1,4,-2],A3=[-3,1,4];b1=[-3,2,4]T,b2=[100,-200,345]T,(2) A行分别为A1=[1,0,8,0.8],A2=[0.8,1,0.8],A3=[0.8,0.8,1];b1=[3,2,1]T,b2=[5,0,-10]T,(3)A行分别为A1=[1,3],A2=[-7,1];b=[4,6]T2.1计算结果初值均为0矩阵带入(1)A行分别为A1=[6,2,-1],A2=[1,4,-2],A3=[-3,1,4];b1=[-3,2,4]T,b2=[100,-200,345]T2) A行分别为A1=[1,0,8,0.8],A2=[0.8,1,0.8],A3=[0.8,0.8,1];b1=[3,2,1]T,b2=[5,0,-10]TT2.2结果分析ρ小于1,故方程组雅可比迭代收第一小题的经计算谱半径为5427B(=).0敛。
计算方法上机实习题
数值计算方法上机实习题1. 设⎰+=105dx xx I nn , (1) 由递推公式nI I n n 151+-=-,从0=0.1822I , 0=0.1823I 出发,计算20I ; (2) 20=0I ,20=10000I , 用nI I n n 51511+-=-,计算0I ;(3) 分析结果的可靠性及产生此现象的原因(重点分析原因)。
解:(1)程序如下: clear all clc I=0.1822; %题中的已知数据 for n=1:20; I=(-5)*I+1/n; %由递推公式所得 end fprintf('I20=%f\n',I) M=0.1823; %与I 的计算结果形成对比for i=1:20; M=(-5)*M+1/i; %由递推公式所得 end fprintf('M20=%f\n',M) 输出结果为: I20=-11592559237.912731 M20=-2055816073.851284 (2)程序如下: clear all clc I=0; %赋予I20的初始值 for n=0:19; I=(-1/5)*I+1/(5*(20-n)); %有递推公式得 end fprintf('I0=%f\n',I)M=10000; for i=0:19; M=(-1/5)*M+1/(5*(20-i));%有递推公式得 end fprintf('M0=%f\n',M) 输出结果为: I0=0.182322 M0=0.182322(3)由输出结果可看出第一种算法为不稳定算法,第二中算法为稳定算法。
由于误差*000***21111120115(5)5()555nn n n n n n n n n e I I e I I I I I I e e e n n------=-=-=-+--+=-===第一种算法为正向迭代算法,每计算一步误差增长5倍,虽然所给的初始值很接近,随着n 的增大,误差也越来越大。
数值计算方法上机实验题
数值计算方法上机实验题数值计算方法上机实验实验内容:1.要求:分别用复化梯形法,复化Simpson 法和 Romberg 公式计算.2.给定积分dx e x31和dx x ?311 ,分别用下列方法计算积分值要求准确到510- ,并比较分析计算时间. 1)变步长梯形法; 2)变步长 Simpson 法; 3) Romberg 方法.算法描述:1、复合梯形法:?=tdt t a t V 0)()( ))()(2)((211∑-=++=n k k n b f x f a f hT输入被积函数数据点t,a. 输出积分值.n T复合Simpson 法:?=tdt t a t V 0)()( ))()(2)(4)((6101121∑∑---=++++=n k n k k k n b f x f x f a f h输入被积函数f(x),积分区间[a,b]和n 输出复合Simpson 积分值n S步1 .);()(;a x b f a f S nab h n ?-?-? 步2 对n k ,,2,1 =执行).(2;2);(4;2x f S S hx x x f S S h x x n n n n +?+?+?+?步3 n n S hS ??6步4 输出n SRomberg 积分法:根据已知数据对其进行多项式拟合得出p(x);f(x)?p(x); 输入被积函数f(x),积分区间端点a,b,允许误差ε 输出 Romberg 积分值n R 2 步1 .0;0;0;0));()((2;1111?===+?-?k R C S b f a f hT a b h 步2 反复执行步3→步9. 步3 .2;0h a x S +步4 反复执行步5→步6. 步5 ;);(h x x x f S S +?+?步6 若x ≥b,则退出本层循环. 步7 执行.6316364;1511516;3134;2212212212212C C R S S C T T S S h T T -?-?-?+?步8 执行.1;;;;;2;2121212112+-?k k R R C C S S T T hh R R e 步9 若e ≤ε且k ≥5,则退出循环. 步10 .22R R n ? 步11 输出.2n R2、变步长梯形算法:功能求积分ba)(dx x f ,允许误差为ε。
数值分析计算实习第一题
直接用定义: ������������(������������)2 = ‖������������‖2‖������������−1‖2
求 A 的条件数很繁琐,需要先进行化简:
首先:
由于 A 是对称矩阵,
‖������������‖2 = �������������max(������������������������������������)
说明 :
1. 在所用的算法中,凡是要给出精度水平的ε,都取 ������������=10−12。
2. 选择算法的时候应使矩阵 A 的所有零元素都不存储。
3. 打印以下内容:
(1)算法设计方案和思路。
(2)全部源程序。
(3)特征值������������1,������������501,������������������������,������������������������������������(������������=1,2,⋯,39)以及������������������������������������������������(������������)2, det������������的值(采用 e 型输出实型数,并 至少显示 12 位有效数字)。
λi[16] -2.533970311130E+00 λi[38] 8.648666065193E+00
λi[17] -2.003230769563E+00 λi[39] 9.254200344575E+00
λi[18] -1.503557611227E+00 cond(A)2 1.925204273903E+03
λi[19] -9.935586060080E-01 det(A) 2.772786141752E+118
数值计算方法第一次上机实习报告
数值计算方法第一次实习报告一、 实习题目。
1. 分别用高斯-赛德尔迭代法、雅克比迭代法、列主元消去法解下列方程组。
⎪⎪⎩⎪⎪⎨⎧=+++=+++=+++=+++ 1.84711.2671x 0.2568x 0.2471x 0.2368x1.74710.2271x 1.2168x 0.2071x 0.1968x 1.64710.1871x 0.1768x 1.1675x 0.1582x 1.54710.1490x 0.1397x 0.1254x 1.1161x4321432143214321 2. 用迭代法求x 5-x-0.2=0的正根,要求精确到小数点后第五位。
二、算法原理。
1、雅可比迭代法基本原理 将矩阵分解为,其中则式可记为,变形可得,可逆时,有于是得到迭代的过程为式中,,即2、高斯-赛德尔迭代法基本原理赛德尔迭代法是对雅可比迭代法的一种改进,雅可比迭代法是在每一步计算的各个分量时均只用到中的分量。
实际上,在计算时,分量都已经计算出来而没有被直接利用,因此可以考虑以来代替计算。
即121311121232222121-1,212121000, ,0000n n n n nn a a a a a a a a D L U aa a a a a a ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥==-=-⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎢⎥⎢⎥⎣⎦⎣⎦=Ax b ()=D L U x b--()=+Dx L U x b+D ()11220nn a a a ≠ ()11=D +x L U x D b --+=+x Bx f ()11D,=D B L U f b --=+()=11=+=1,2,,n i i ij j j ii j i x b a x i n a ≠⎛⎫ ⎪ ⎪ ⎪⎝⎭∑ ()+1k x ()k x ()+1k i x ()()+1+1-1,,k k i i x x ()()+1+1-1,,k k i i x x ()()1-1,,kki x x矩阵形式为,可得,于是赛德尔迭代法的矩阵形式为式中,。
数值计算方法上机实习题
数值计算方法上机实习题(总13页)--本页仅作为文档封面,使用时请直接删除即可----内页可以根据需求调整合适字体及大小--数值计算方法上机实习题1.设⎰+=105dx xx I nn , (1) 由递推公式nI I n n 151+-=-,从I 0=0.1824, 0=0.1823I 出发,计算20I ; (2) 20=0I ,20=10000I , 用nI I n n 515111+-=--,计算0I ;(3) 分析结果的可靠性及产生此现象的原因(重点分析原因)。
答:第一个算法可得出e 0=|I 0−I 0∗| e n =|I n −I n ∗|=5n |e 0|易知第一个算法每一步计算都把误差放大了5倍,n 次计算后更是放大了5n 倍,可靠性低。
第二个算法可得出e n =|I n −I n ∗| e 0=(1)n|e n |可以看出第二个算法每一步计算就把误差缩小5倍,n 次后缩小了5n 倍,可靠性高。
2.求方程0210=-+x e x 的近似根,要求41105-+⨯<-k k x x ,并比较计算量。
(1) 在[0,1]上用二分法; (1) [0,1]上的二分法二分法子程序:function [root,n]=EFF3(f,x1,x2) %第二题(1)二分法f1=subs(f,symvar(f),x1);%函数在x=x1的值 f2=subs(f,symvar(f),x2);%x=x2 n=0;%步数er=5*10^-4;%误差 if(f1==0) root=x1; return; elseif(f2==0) root=x2; return;elseif(f1*f2>0)disp('两端点函数值乘积大于0!'); return; elsewhile(abs(x1-x2)>er)%循环 x3=(x1+x2)/2;f3=subs(f,symvar(f),x3); n=n+1; if(f3==0) root=x3; break;elseif(f1*f3>0) x1=x3; elsex2=x3; end endroot=(x1+x2)/2;%while 循环少一步需加上 end计算根与步数程序:fplot(@(x) exp(x)+10*x-2,[0,1]); grid on; syms x;f=exp(x)+10*x-2; [root,n]=EFF3(f,0,1);fprintf('root=% ,n=%d \n',root,n);计算结果显示:root= ,n=11(2) 取初值00=x ,并用迭代1021x k e x -=+;(2) 初值x 0=0迭代 迭代法子程序:function [root,n]=DDF(g,x0,err,max) (接下页) %root 根,n+1步数,g 函数,x0初值,err 误差,max 最大迭代次数 X(1)=x0; for n=2:maxX(n)=subs(g,symvar(g),X(n-1)); c=abs(X(n)-X(n-1)); root=X(n); if(c<err)计算根与步数程序:syms x;f=(2-exp(x))/10; (接下页) x0=0;err=5*10^(-4); max=100;[root,n]=DDF(f,x0,err,max); fprintf('root=% ,n=%d \n',root,n);计算结果显示: root= ,n=4(3) 加速迭代的结果;(4) 取初值00 x ,并用牛顿迭代法;(5) 分析绝对误差。
[整理]-8-《数值计算方法》实习作业(模板-小)
2.1函数图形与极限2.1.1 实验目的1.熟悉Mathematica 基本绘图语句。
2.掌握函数极限的有关操作命令。
3.学会利用Mathematica 软件对函数进行分析研究。
4.熟悉Mathematica 二元函数绘图语句。
2.1.2 实验内容【基本语句】1.Plot[f[x],{x,xmin,xmax},选项]; 功能: 画出函数f[x] 从min 到max 间的图形;2.Plot[{f1[x],f2[x],...},{x,xmin,xmax},选项]; 功能: 在同一坐标系下画出函数f1,f2,...的图形。
3. ParametricPlot[{fx,fy},{t,tmin,tmax}]; 功能: 画出参数方程fx=x(t),fy=y(t)的图形;ParametricPlot[{{f1x,f1y},{f2x,f2y}},{t,tmin,tmax}]; 功能:在同一坐标系下画出用参数方程表示的两幅函数图形。
【备注】fx,fy 的给出方式:⑴fx=x(t) , fy=y(t)⑵fx=x ,fy=f(x)与fx=f(x) ,fy=x 构成反函数的图形关系⑶r=r(t) , fx=r(t)Cos(t) , fy=r(t)Sin(t)4. Show[tu1,tu2] 功能:将tu1及tu2两幅函数图形重叠在一起,将两个函数图形一起显示。
5. Plot3D[f[x,y],{x,x0,x1},{y,y0,y1}] 功能:作出函数f[x,y]在区域[x0,x1]×[y0,y1]上的图形; ParametricPlot3D[{x[u,v],y[u,v],z[u,v]},{u,u0,u1},{v,v0,v1}] 功能:作出参数方程表示的曲面。
6. Limit[f[x],x->x0] 功能:求函数f[x]在x0处的极限。
7. Limit[f[x],x->x0,Direction->+1] 功能:求函数f[x]在x0处的左极限。
数值计算方法上机实习题答案.doc
1.设I n 1 x ndx ,0 5 x( 1)由递推公式 I n 5I n 11,从 I 0的几个近似值出发,计算I 20;n解:易得: I 0 ln6-ln5=0.1823, 程序为:I=0.182;for n=1:20I=(-5)*I+1/n;endI输出结果为: I 20= -3.0666e+010( 2)粗糙估计 I 20,用 I n 1 1I n 1 1 ,计算 I 0;5 5n0.0079 1 x 20 1 x 200.0095因为dx I 20dx 6 5所以取 I 20 1(0.0079 0.0095) 0.0087 2程序为: I=0.0087;for n=1:20I=(-1/5)*I+1/(5*n);endII 0= 0.0083( 3)分析结果的可靠性及产生此现象的原因(重点分析原因 )。
首先分析两种递推式的误差;设第一递推式中开始时的误差为E0 I 0 I 0,递推过程的舍入误差不计。
并记 E n I n I n,则有 E n 5E n 1 ( 5) n E0。
因为 E20( 5) 20 E0 I 20,所此递推式不可靠。
而在第二种递推式中E0 1E1 (1)n E n,误差在缩小,5 5所以此递推式是可靠的。
出现以上运行结果的主要原因是在构造递推式过程中,考虑误差是否得到控制,即算法是否数值稳定。
2.求方程e x10x 2 0 的近似根,要求x k 1x k 5 10 4,并比较计算量。
(1)在 [0, 1]上用二分法;程序: a=0;b=1.0;while abs(b-a)>5*1e-4c=(b+a)/2;if exp(c)+10*c-2>0b=c;else a=c;endendc结果: c =0.0903( 2)取初值x0 0,并用迭代 x k 1 2 e x ;10程序: x=0;a=1;while abs(x-a)>5*1e-4a=x;x=(2-exp(x))/10;endx结果: x =0.0905(3)加速迭代的结果;程序: x=0;a=0;b=1;while abs(b-a)>5*1e-4a=x;y=exp(x)+10*x-2;z=exp(y)+10*y-2;x=x-(y-x)^2/(z-2*y+x);b=x;endx结果: x =0.0995( 4)取初值x00 ,并用牛顿迭代法;程序: x=0;a=0;b=1;while abs(b-a)>5*1e-4a=x;x=x-(exp(x)+10*x-2)/(exp(x)+10); b=x;end x 结果: x =0.0905( 5) 分析绝对误差。
数值计算方法上机实习题NEW
数值计算方法上机实习题1. 设,(1) 由递推公式,从 , 出发,计算 ; 程序如下:function I=myhs(I0,n) if n>=1I=myhs(I0,n-1)*(-5)+1/n; elseif n==0 I=I0; end命令行窗口输入: I0=0.1822;I1=myhs(I0,20); I0=0.1823;I2=myhs(I0,20);运行结果:当I0=0.1822时,I20=-1.1593e+10。
当I0=0.1823时,I20= -2.0558e+9。
(2) ,20=10000I , 用,计算0I ; 程序如下:function I=myhs2(I20,n) if (n<20)I=myhs2(I20,n+1)*(-1)/5+1/(5*(n+1)); elseif n==20 I=I20; end命令行窗口输入:I20=0;I1=myhs2(I20,20); I20=10000;I2=myhs2(I20,20);运行结果:当I 20=0时,I 0=0.182321556793955。
当I 20=10000时,I 0= 0.182321556898812。
(3) 分析结果的可靠性及产生此现象的原因(重点分析原因)。
根据上式,假设I n 的真值为I*,误差为e n ,即e=I*-I n 。
综合递推式,有e n =-5*e n-1,这意味着哪怕开始只有一点点误差,只要n 足够大,按照这种计算一步误差增长5倍的方式,所得的结果总是不可信的,因此整个算法是不稳定的。
根据上式,假设I n 的真值为I*,误差为e n ,即e=I*-I n 。
综合递推式,有e n-1. =(-1/5)*e n ,按照这种计算误差会以每步缩小到1/5的方式进行,根据(2)得到的结果而言,该算法是相对稳定的。
2. 求方程0210=-+x e x 的近似根,要求41105-+⨯<-k k x x ,并比较计算量。
数值方法计算实习题1
信计091 龚立丽200900901004数值方法计算实习题要求:1、用Matlab语言或你熟悉的其他算法语言编写程序,使之尽可能具有通用性;2、根据上机计算实践,对所使用的数值方法的特点、性质、有效性、误差和收敛性等方面进行必要的讨论和分析;3、完成计算后写出实验报告,内容包括:课题名称、解决的问题、采用的数值方法、算法程序、数值结果、对实验结果的讨论和分析等;4、特别说明:严禁抄袭,否则一经发现,所有雷同实验报告最多评为及格。
一、下表给出了飞行中鸭子的上部形状的节点数据,试用三次样条插值函数(自然边界条件)和20次Lagrange插值多项式对数据进行插值。
用图示出给定的数据,以及()s x和20()L x。
x0.9 1.3 1.9 2.1 2.6 3.0 3.9 4.4 4.7 5.0 6.0y 1.3 1.5 1.85 2.1 2.6 2.7 2.4 2.15 2.05 2.1 2.25x7.0 8.0 9.2 10.5 11.3 11.6 12 12.6 13.0 13.3y 2.3 2.25 1.95 1.4 0.9 0.7 0.6 0.5 0.4 0.25解:>> x=[0.9 1.3 1.9 2.1 2.6 3.0 3.9 4.4 4.7 5.0 6.0 7.0 8.0 9.2 10.5 11.3 11.6 12 12.6 13.0 13.3];>> y=[1.3 1.5 1.85 2.1 2.6 2.7 2.4 2.15 2.05 2.1 2.25 2.3 2.25 1.95 1.4 0.9 0.7 0.6 0.5 0.4 0.25];(1)三次样条插值法在MATLAB中编写m文件function[f,f0]=scyt(x,y,y2_1,y2_N,x0)%y2_1和y2_N分别为自然边界条件%插值点x的坐标:x0syms t;f=0.0;f0=0.0;if(length(x)==length(y))n=length(x);elsedisp('x和y的维数不相等');return;endfor i=1:nif(x(i)<=x0)&&(x(i+1)>=x0)index=i;break;endendA=diag(2*ones(1,n));A(1,2)=1;A(n,n-1)=1;u=zeros(n-2,1);lamda=zeros(n-1,1);c=zeros(n,1);for i=2:n-1u(i-1)=(x(i)-x(i-1))/(x(i+1)-x(i-1));lamda(i)=(x(i+1)-x(i))/(x(i+1)-x(i-1));c(i)=3*lamda(i)*(y(i)-y(i-1))/(x(i)-x(i-1))+3*u(i-1)*(y(i+1)-y(i))/(x(i+1)-x(i) );A(i,i+1)=u(i-1);A(i,i-1)=lamda(i);endc(1)=3*(y(2)-y(1))/(x(2)-x(1))-(x(2)-x(1))*y2_1/2;c(n)=3*(y(n)-y(n-1))/(x(n)-x(n-1))-(x(n)-x(n-1))*y2_N/2;m=zgf(A,c);h=x(index+1)-x(index);f=y(index)*(2*(t-x(index))+h)*(t-x(index+1))^2/h/h/h+y(index+1)*(2*(x(index+1)-t)+h)*(t-x(index))^2/h/h/h+m(index)*(t-x(index))*(x(index+1)-t)^2/h/h-m(index+1 )*(x(index+1)-t)*(t-x(index))^2/h/h;f0=subs(f,'t',x0);其中的zgf函数为追赶法,其程序为function x=zgf(A,b)n = rank(A);for(i=1:n)if(A(i,i)==0)disp('Error: 对角有元素为0!');return;endend;d = ones(n,1);a = ones(n-1,1);c = ones(n-1);for(i=1:n-1)a(i,1)=A(i+1,i);c(i,1)=A(i,i+1);d(i,1)=A(i,i);endd(n,1) = A(n,n);for(i=2:n)d(i,1)=d(i,1) - (a(i-1,1)/d(i-1,1))*c(i-1,1);b(i,1)=b(i,1) - (a(i-1,1)/d(i-1,1))*b(i-1,1);endx(n,1) = b(n,1)/d(n,1);for(i=(n-1):-1:1)x(i,1) = (b(i,1)-c(i,1)*x(i+1,1))/d(i,1); end在MATLAB 中输入指令 >> [f,f0]=scyt(x,y,0,0) f =1000/729*(27/5*t-1377/100)*(t-39/10)^2+1000/729*(522/25-24/5*t)*(t-3)^2+100/81*(-6396162352027119/288230376151711744*t+19188487056081357/288230376151711744)*(39/10-t)^2-100/81*(-176836856862157557/90071992547409920+4534278381080963/9007199254740992*t)*(t-3)^2 f0 =2.5851 得三次样条插值函数 S(x)=1000/729*(27/5*x-1377/100)*(x-39/10)^2+1000/729*(522/25-24/5*x)*(x-3)^2+100/81*(-6396162352027119/288230376151711744*x+19188487056081357/288230376151711744)*(39/10-x)^2-100/81*(-176836856862157557/90071992547409920+4534278381080963/9007199254740992*x)*(x-3)^2>> xi=0.9:0.01:13.3;yi=interp1(x,y,xi,'spline'); >> title('试验一--三次样条插值图示')024********0.511.522.53试验一--三次样条插值图示(2)用拉格朗日法插值 %定义Lagrange 程序function f=Language(x,y,x0) syms t ;if (length(x)==length(y)) n=length(x);elsedisp('x 和y 的维数不相等'); return ; end f=0.0; for (i=1:n) l=y(i); for (j=1:i-1)l=l*(t-x(j))/(x(i)-x(j)); end ;for (j=i+1:n)l=l*(t-x(j))/(x(i)-x(j)); end ; f=f+l; simplify(f); if (i==n)if (nargin==3) f=subs(f,'t',x0); elsef=collect(f); f=vpa(f,6); end end end>> Language(x,y) ans =52462.6*t+189995.*t^3-189851.*t^4+136778.*t^5-11.3161*t^12-.277283e-6*t^18+1.18284*t^13-73866.6*t^6+.111076e-4*t^17-.976904e-1*t^14+.427949e-8*t^19-.307453e-10*t^20+30677.6*t^7+2564.20*t^9-9968.98*t^8+.628590e-2*t^15-525.813*t^10-9652.78-.308159e-3*t^16+86.2514*t^11-128683.*t^2二、已知Wilson 矩阵1078775658610975910A ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦,且向量32233331b ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦,则方程组A x b =有准确解[]1111Tx =。
——数值分析上机题
.......................课程名称:数值分析上机实习报告姓名:学号:专业:联系电话:目录序言 (3)第1章必做题 (4)1.1必做题第一题 (4)1.1.1题目 (4)1.1.2 分析 (4)1.1.3 计算结果 (4)1.1.3 总结 (6)1.2必做题第二题 (6)1.2.1题目 (6)1.2.2分析 (6)1.2.3计算结果 (6)1.2.4结论 (8)1.1必做题第一题....................................................................... 错误!未定义书签。
1.1.1题目 ............................................................................ 错误!未定义书签。
第2章选做题 (8)2.1选做题第一题 (8)2.1.1题目 (8)2.1.2分析 (8)2.1.3计算结果 (8)附录 (10)附录一:必做题第一题程序 (10)附录二:必做题第二题程序 (11)附录三:选做题第一题的程序 (13)序言本次数值分析上机实习采用Matlab数学软件。
Matlab是一种用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境。
在数值分析应用中可以直接调用Matlab软件中已有的函数,同时用户也可以将自己编写的实用程序导入到Matlab函数库中方便自己调用。
基于Matlab数学软件的各种实用性功能与优点,本次数值分析实习决定采用其作为分析计算工具。
1.编程效率高MATLAB是一种面向科学与工程计算的高级语言,允许使用数学形式的语言编写程序,且比BASIC、FORTRAN和C等语言更加接近我们书写计算公式的思维方式,用MATLAB编写程序犹如在演算纸上排列出公式与求解问题。
因此,MATLAB语言也可通俗地称为演算纸式科学算法语言。
数值计算方法上机实习题考证
---------------------------------------------------此文档包含我们计算方法的经典算法包含(数值计算方法上机实习题)1. 设⎰+=105dx x x I nn , (1) 由递推公式nI I n n 151+-=-,从0I 的几个近似值出发,计算20I ; (2) 粗糙估计20I ,用n I I n n 51511+-=-,计算0I ; (3) 分析结果的可靠性及产生此现象的原因(重点分析原因)。
(1) 解答:n=0,0.1823)05ln()15ln()5(51515101010=+-+=++=+=+=⎰⎰⎰x d x dx x dx x x I nn 这里可以用for 循环,while 循环,根据个人喜好与习惯:for 循环程序: While 循环程序:I=0.1823; I=0.1823;for n=1:20 i=1;I=(-5)*I+1/n; while i<21End I=(-5)*I+1/i;I i=i+1;fprintf('I20=%f',I) endI = -2.0558e+009 >> II20=-2055816073.851284>> I = -2.0558e+009(2) 粗略估计I 20: Mathcad 计算结果: for 循环程序: While 循环程序: >> I=0.007998; I=0.007998;>> for n=1:20 n=1;I=(-0.2)*I+1/(5*n); while n<21End I=(-0.2)*I+1/(5*n);>> I n=n+1;I =0.0083 end>> II =0.0083(3) 算法误差分析:计算在递推过程中传递截断误差和舍入误差第一种算法:(从1——>20)01x x 205x +⎛⎜⎜⎜⎠d 7.998103-⨯=*000e I I =-***21111120115(5)5()555n n n n n n n n n n e I I I I I I e e e n n------=-=-+--+=-===误差放大了5n 倍,算法稳定性很不好;第二种算法:(从20——>1)*n n ne I I =- ***111111111()()555555n n n n n n n n e I I I I I I e n n ---=-=-+--+=-=0111...()55n n e e e === 误差在逐步缩小,算法趋近稳定,收敛。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
数值计算方法上机实习题1 1(2)丨20=0, 120 = 10000,用 I n 」I n 」 一,计算丨0 ; 5 5n答:第一个算法可得出易知第一个算法每一步计算都把误差放大了5倍,n 次计算后更是放大了 5n倍,可靠性低。
第二个算法可得出(1 )由递推公式I n = -51心-,从nI o =O.1823出发,计算* ;(3)分析结果的可靠性及产生此现象的原因重点分析原因。
可以看出第二个算法每一步计算就把误差缩小5倍,n次后缩小了5n倍,可靠性高。
2.求方程e^+IOx —2=0的近似根,要求 x k *—x k <^10-,并比较计算量。
(1)在[0 , 1]上用二分法;(1) [0, 1]上的二分法二分法子程序:function [root,n]=EFF3(f,x1,x2) %第二题⑴二分法f1= subs(f,symvar(f),x1);% 函数在 x=x1 的值 f2=subs(f,symvar(f),x2);%x=x2 n=0;%步数 er=5*10A-4;% 误差 if(f1==0)root=x1; return;elseif(f2==0)root=x2; return;elseif(f1*f2>0)disp('两端点函数值乘积大于 0!');return;elsewhile(abs(x1-x2)>er)% 循环 x3=(x1+x2)/2;f3=subs(f,symvar(f),x3);n=n+1; if(f3==0)root=x3; break; elseif(f1*f3>0)x1=x3; elsex2=x3; end endroot=(x1+x2)/2;%while 循环少一步需加上 end、2 _e (2) 取初值x 0 =0,并用迭代x k d :10(2)初值x o =O 迭代计算根与步数程序:迭代法子程序:syms x;function [root,n]=DDF(g,x0,err,max) (接下页) f=(2-exp(x))/10;(接下页)计算根与步数程序:fplot(@(x) exp(x)+10*x-2,[0,1]); grid on; syms x; f=exp(x)+10*x-2; [root,n]=EFF3(f,0,1);fprintf('root=%6.8f ,n=%d \n',root,n);计算结果显示:root=0.09057617 ,n=11(3)加速迭代的结果(4)取初值X。
=0,并用牛顿迭代法;(5)分析绝对误差。
答:可以看到,在同一精度下,二分法运算了11次,题设迭代算式下运算了4次,加速迭代下运算了2次,牛顿迭代下运算了2次。
因不动点迭代法和二分法都是线性收敛的,但二分法压缩系数比题设迭代方法大,收敛速度较慢。
加速迭代速度是超线性收敛,牛顿法是二阶,收敛速度快。
3.钢水包使用次数多以后,钢包的容积增大,数据如下:ax + b试从中找出使用次数和容积之间的关系,计算均方差。
(用y 拟合)c + x拟合曲线程序:x=2:16;y=[6.42 8.2 9.58 9.5 9.7 10 9.93 9.99 10.49 10.5910.60 10.8 10.6 10.9 10.76];[f,fval]=fminsearch(@delta1,[0,0,0],optimset,x,y)J%fminsearch为求解多元函数的最小值函数%f为多元函数初值x0附近的极小值点%fval为极小值k=2:100;y1= (f(1).*k+f(2))./(f(3)+k); figure1=figure('Color',[1 1 1]);gxt=plot(x,y,'xr',k,y1,'k-');legend(原数据’,'拟合曲线','Location','northwest'); %y 为数据点连接曲线,y1为拟合曲线title('函数y=(ax+b)/(x+c)的拟合','FontSize',14,'FontWeight','Bold'); xlabel('次数','FontSize',14,'FontWeight','Bold'); ylabel('容积','FontSize',14,'FontWeight','Bold');set(gxt,'LineWidth',1.5);grid on;%计算均方差for i=1:15 构造函数子程序:function delta=delta1(f,x,y)a=f(1);b=f(2);c=f(3);delta=0;for k=1:15delta=delta+((x(k)+c)*y(k)-(a*x(k)+b))A2;end计算结果显示:拟合出的方程为:(x+-0.7110)y=11.2657x+-15.5024均方差为033165089总结:指标选择,因题设方程I & E需[片—签普「为非线性的,要转化为线性方程故需提指标为:y2(i)=(f(1).*x(i)+f(2))./(f(3)+x(i)); (接下页)end j=0; for i=1:15j=i+(y(i)-y2(i))A2; endjfc=sqrt(j/15);fprintf('拟合出的方程为:(x+%4.4f)y=%4.4fx+%4.4f \n 均方差 为:%4.8f \n',f(3),f(1),f(2),jfc);i=i亦2+ l)y£— oxf - 力]二1=1(cx^-i)y=ax+b f l= J-jJfcXj +1)给一叫 一 b] 其驻点方程为:、亠£ =》2[宓+Dyj-oXi- >](-%!)= 0;i=i计算结果显示:12610980 10 20 30 40 50 60 70 80 90 100次数『4-1 0 -1 0 0、『0、-1 4 -1 0 -1 050 -1 4 -1 0 -1 ,b = -2 -1 0 -1 4 -1 050 -1 0 -1 4 -1-2 1° 0 -1 0 -1 4」3」(1) JACOBI 迭代;(2) GAUSS-SEIDEL 迭代;分析下列迭代法的收敛性,并求Xk 1 -xk<10^的近似解及相应的迭代次数。
(3)SOR迭代(取=0.1:0.1:1.9,找到迭代步数最少的* )。
逆幕迭代法子程序:function [lamda,V ,k]=NMSF2(A,v,eps,p) [n,n]=size(A); A=A-p*eye( n); v0=v;[tmax,tindex]=max(abs(v0)); lamd0=v0(tindex); u0=v0/lamd0; flag=0; k=0;while(flag==0)V=A\u0;[tmax,tindex]=max(abs(V)); lamd1=V(tindex); u0=V/lamd1;if (abs((lamd0)A(-1)-(lamd1)A(-1)))v=eps flag=1; endlamd0=lamd1; k=k+1;end (接下页) lamda=(lamd1)A(-1)+p; V=u0';逆幕迭代计算程序:A=[6 3 1;3 2 1;1 1 1]; v=[1 1 1]'; eps=0.001; p=11;[lamda,U,k]=NMSF2(A,v,eps,p);fprintf('最接近11的特征值为:%4.8f\n 特征向量: \n%4.8f\n %4.8f\n%4.8f\n',lamda,U(1),U(2),U(3));计算结果显示:最接近于11的特征值为:7.87313712 特征向量:1.00000000 0.54922698 0.225456186.用经典R-K 方法求解初值问题"6 3r5•用逆幕迭代法求 A =3 2 1 最接近于1 1 111的特征值和特征向量,准确到10’。
'y; = _2y<i + y2 +2sin x”2 = yi —2y2 +2cosx—2sin x和精确解丿y i(x) =2e」Minx比较,进行误差分析得到结论,图形显示精确解和数值解。
y2(x) = 2e» +cosx经典四阶R-K方法R-K程序:f=@(x,y)[-2*y(1)+y(2)+2*sin(x);y(1)-2*y(2)+2*cos(x)-2*sin(x)];g=@(x,y)【-2*y(1)+y(2)+2*sin(x);998*y(1)-999*y(2)+999*cos(x)-999*sin(x)]; (接下页)x=0:0.1:10;Y1= 2*exp(-x)+sin(x);Y2=2*exp(-x)+cos(x);x仁0:10/10000:10;Y3=2*exp(-x1)+sin(x1);Y4=2*exp(-x1)+cos(x1);N=100;N 1=10000;[a1,b1]=ode23(f,[0:0.1:10],[2 3]);[a2,b2]=ode23(g,[0:1/1000:10],[2 3]);[a3,b3]=ode45(f,[0:0.1:10],[2 3]);[a4,b4]=ode45(g,[0:1/1000:10],[2 3]);n=length(b1);figure1=figure('Color',[1 1 1]);dbt=plot(x,b1(:,1),'r-',x,b1(:,2),'k-',x,Y1,'g--',x,Y2,'y--'); legend('y1','y2','精确解Y1','精确解Y2');title('四阶RK算法下的方程组的数值解&方程的精确解','FontSize',14,'FontWeight','Bold');ylabel( Y,'FontSize',14,'FontWeight','Bold');xlabel('X','FontSize',14,'FontWeight','Bold');set(dbt,'LineWidth',1.5);grid on;Q=b1(:,1)-Y1';W=b1(:,2)-Y2';Q1= b2(:,1)-Y3';W 仁b2(:,2)-Y4';Q2=b3(:,1)-Y1';W2=b3(:,2)-Y2';Q3=b4(:,1)-Y3';W3=b4(:,2)-Y4';ER1=max(abs(Q));ER2=max(abs(W));ER3=max(abs(Q1));ER4=max(abs(W1));fprintf('\n (1)中y1在ode23下求出的最大误差为:%4.8f\n( 1 )中y2在ode23下求出的最大误差为:%4.8f\n',ER1,ER2);fprintf('\n (2)中y1在ode23下求出的最大误差为:%4.8f\n(2)中y2在ode23下求出的最大误差为:%4.8f\n',ER3,ER4);fprintf('\n (1)中y1在ode45下求出的最大误差为:%4.8f\n ( 1 )中y2在ode45下求出的最大误差为:%4.8f\n',ER5,ER6);fprintf('\n (2)中y1在ode45下求出的最大误差为:%4.8f\n(2)中y2在ode45下求出的最大误差为:%4.8f\n',ER7,ER8);计算与显示程序:(1 )中y1在ode23下求出的最大误差为:0.00169850(1 )中y2在ode23下求出的最大误差为:0.00207422(2)中y1在ode23下求出的最大误差为:0.00000583(2)中y2在ode23下求出的最大误差为:0.00581329(1 )中y1在ode45下求出的最大误差为:0.00052155(1 )中y2在ode45下求出的最大误差为:0.00045793(2)中y1在ode45下求出的最大误差为:0.00000276(2)中y2在ode45下求出的最大误差为:0.00275331总结:对比2种方程在ode23,ode45命令下,与精确值的误差可以发现,在此题中,ode45的误差相对于ode23的误差小,较适合用在此类钢性问题中。