数值分析(王能超)matlab上机作业
数值分析上机作业(MATLAB)
将系数矩阵 A 分解为：A=L+U+D
Ax=b
⇔ (D + L +U)x = b ⇔ Dx = −(L + U )x + b ⇔ x = −D −1(L + U )x + D −1b x(k +1) = −D −1 (L + U ) x(k ) + D −1b
输入 A，b 和初始向量 x
迭代矩阵 BJ , BG
否
ρ(B) < 1?
按雅各比方法进行迭代
否
|| x (k+1) − x(k) ||< ε ?
按高斯-塞德尔法进行迭代
否
|| x(k+1) − x (k ) ||< ε ?
输出迭代结果
图 1 雅各布和高斯-赛德尔算法程序流程图
1.2 问题求解
按图 1 所示的程序流程，用 MATLAB 编写程序代码，具体见附录 1。解上述三个问题 如下
16
-0.72723528355328
0.80813484897616
0.25249261987171
17
-0.72729617968010
0.80805513082418
0.25253982509100
18
-0.72726173942623
0.80809395746552
0.25251408253388
0.80756312717373
8
-0.72715363032573
0.80789064377799
9
-0.72718652854079
数值分析上机题3
数值分析上机题目3实验一1．根据Matlab 语言特点，描述Jacobi 迭代法、Gauss-Seidel 迭代法和SOR 迭代法。
2．编写Jacobi 迭代法、Gauss-Seidel 迭代法和SOR 迭代法的M 文件。
3．给定2020⨯∈R A 为五对角矩阵⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡------------------321412132141412132141412132141412132141213O O O O (1)选取不同的初始向量)0(x 及右端面项向量b ，给定迭代误差要求，分别用编写的Jacobi 迭代法和Gauss-Seidel 迭代法程序求解，观察得到的序列是否收敛？若收敛，通过迭代次数分析计算结果并得出你的结论。
(2)用编写的SOR 迭代法程序， 对于(1)所选取的初始向量)0(x 及右端面项向量b 进行求解，松驰系数ω取1<ω<2的不同值，在5)1()(10-+≤-k k x x 时停止迭代，通过迭代次数分析计算结果并得出你的结论。
实验11、 根据MATLAB 语言特点，描述Jacobi 迭代法，Gauss-Seidel 迭代法和SOR 迭代法。
2、 编写Jacobi 迭代法，Gauss-Seidel 迭代法和SOR 迭代法的M 文件。
Jacobi 迭代法function [x1,k]=GS_2(A,b)n=length(A);D=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);x1=zeros(n,1);x0=3*ones(n,1);k=0; while norm(x1-x0,1)>10^(-7)&k<100 k=k+1;x0=x1;x1=D\((L+U)*x0+b);endk=kx=x1Gauss-Seidel迭代法function [x1,k]=GS_h(A,b)n=length(A);D=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);x1=zeros(n,1);x0=3*ones(n,1);k=0; while norm(x1-x0,1)>10^(-7)&k<100 k=k+1;x0=x1;x1=(D-L)\U*x0-D\b;endk=kx=x1SOR迭代法function [x1,k]=SOR_h(A,b)n=length(A);D=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);x1=zeros(n,1);x0=3*ones(n,1);k=0;w=0.96;while norm(x1-x0,1)>10^(-7)&k<100k=k+1;x0=x1;x1=(D-w*U)\(((1-w)*D+w*L)*x0+w*b);endk=kx=x13、采用Jacobi迭代法，Gauss-Seidel迭代法求解五对角矩阵clear,clcA=diag(3*ones(20,1))+diag((-0.5)*ones(19,1),-1)+diag((-0.5)*ones(19,1 ),1)+diag((-0.25)*ones(18,1),-2)+diag((-0.25)*ones(18,1),2);b=sum(A')';[x1,k1]=Jacob_h(A,b)[x2,k2]=GS_h(A,b)运行结果：两种方法都收敛，k1=27,k2=13。
数值计算上机实习题目（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)；给定不同的初始值，观察用牛顿法求解时所需要的迭代步数，并与上面第一题的迭代步数比较。
数值分析上机实践报告
数值分析上机实践报告一、实验目的本次实验主要目的是通过上机操作，加深对数值分析算法的理解，并熟悉使用Matlab进行数值计算的基本方法。
在具体实验中，我们将实现三种常见的数值分析算法：二分法、牛顿法和追赶法，分别应用于解决非线性方程、方程组和线性方程组的求解问题。
二、实验原理与方法1.二分法二分法是一种常见的求解非线性方程的数值方法。
根据函数在给定区间端点处的函数值的符号，不断缩小区间的长度，直到满足精度要求。
2.牛顿法牛顿法是求解方程的一种迭代方法，通过构造方程的泰勒展开式进行近似求解。
根据泰勒展式可以得到迭代公式，利用迭代公式不断逼近方程的解。
3.追赶法追赶法是用于求解三对角线性方程组的一种直接求解方法。
通过构造追赶矩阵，采用较为简便的向前追赶和向后追赶的方法进行计算。
本次实验中，我们选择了一组非线性方程、方程组和线性方程组进行求解。
具体的实验步骤如下：1.调用二分法函数，通过输入给定区间的上下界、截止误差和最大迭代次数，得到非线性方程的数值解。
2.调用牛顿法函数，通过输入初始迭代点、截止误差和最大迭代次数，得到方程组的数值解。
3.调用追赶法函数，通过输入追赶矩阵的三个向量与结果向量，得到线性方程组的数值解。
三、实验结果与分析在进行实验过程中，我们分别给定了不同的参数，通过调用相应的函数得到了实验结果。
下面是实验结果的汇总及分析。
1.非线性方程的数值解我们通过使用二分法对非线性方程进行求解，给定了区间的上下界、截止误差和最大迭代次数。
实验结果显示，根据给定的输入，我们得到了方程的数值解。
通过与解析解进行比较，可以发现二分法得到的数值解与解析解的误差在可接受范围内，说明二分法是有效的。
2.方程组的数值解我们通过使用牛顿法对方程组进行求解，给定了初始迭代点、截止误差和最大迭代次数。
实验结果显示，根据给定的输入，我们得到了方程组的数值解。
与解析解进行比较，同样可以发现牛顿法得到的数值解与解析解的误差在可接受范围内，说明牛顿法是有效的。
数值分析matlab上机实验报告
数值分析matlab上机实验报告matlab软件实验报告数学上机课实验报告matlab实验报告总结数值分析试卷篇一：《MATLAB与数值分析》第一次上机实验报告标准实验报告（实验）课程名称学生姓名：李培睿学号：2013020904026指导教师：程建一、实验名称《MATLAB与数值分析》第一次上机实验二、实验目的1. 熟练掌握矩阵的生成、加、减、乘、除、转置、行列式、逆、范数等运算操作。
（用.m文件和Matlab函数编写一个对给定矩阵进行运算操作的程序）2. 熟练掌握算术符号操作和基本运算操作，包括矩阵合并、向量合并、符号转换、展开符号表达式、符号因式分解、符号表达式的化简、代数方程的符号解析解、特征多项式、函数的反函数、函数计算器、微积分、常微分方程的符号解、符号函数的画图等。
（用.m 文件编写进行符号因式分解和函数求反的程序）3. 掌握Matlab函数的编写规范。
4、掌握Matlab常用的绘图处理操作，包括：基本平面图、图形注释命令、三维曲线和面的填充、三维等高线等。
（用.m 文件编写在一个图形窗口上绘制正弦和余弦函数的图形，并给出充分的图形注释）5. 熟练操作MATLAB软件平台，能利用M文件完成MATLAB的程序设计。
三、实验内容1. 编程实现以下数列的图像，用户能输入不同的初始值以及系数。
并以x，y为坐标显示图像x(n+1) = a*x(n)-b*(y(n)-x(n) ); y(n+1) = b*x(n)+a*(y(n)-x(n) )2. 编程实现奥运5环图，允许用户输入环的直径。
3. 实现对输入任意长度向量元素的冒泡排序的升序排列。
不允许使用sort函数。
四、实验数据及结果分析题目一：①在Editor窗口编写函数代码如下：并将编写的函数文件用“draw.m”储存在指定地址；②在Command窗口输入如下命令：③得到图形结果如下：题目二：①在Editor窗口编写函数代码如下：并将编写的函数文件用“circle.m”储存在指定地址；②再次在Editor窗口编写代码：并将编写的函数文件用“Olympic.m”储存在指定地址；③在Command窗口输入如下指令（半径可任意输入）：④按回车执行，将在图形窗口获得五环旗：题目三：①在Editor窗口编写函数代码如下：并用.将编写的函数文件用“qipaofa.m”储存在指定地址；②在Command窗口输入一组乱序数值，则可以得到升序排序结果如下：五、总结及心得体会1. 要熟悉MATLAB编译软件的使用方法，明白有关语法，语句的基本用法，才可以在编写程序的时候游刃有余，不至于寸步难行。
数值分析作业-matlab上机作业
数值分析———Matlab上机作业学院：班级：老师：姓名：学号：第二章解线性方程组的直接解法第14题【解】1、编写一个追赶法的函数输入a，b，c，d输出结果x，均为数组形式function x=Zhuiganfa(a,b,c,d)%首先说明：追赶法是适用于三对角矩阵的线性方程组求解的方法，并不适用于其他类型矩阵。
%定义三对角矩阵A的各组成单元。
方程为Ax=d%b为A的对角线元素(1~n)，a为-1对角线元素(2~n)，c为+1对角线元素(1~n-1)。
% A=[2 -1 0 0% -1 3 -2 0% 0 -2 4 -3% 0 0 -3 5]% a=[-1 -2 -3];c=[-1 -2 -3];b=[2 3 4 5];d=[6 1 -2 1];n=length(b);u(1)=b(1);y(1)=d(1);for i=2:nl(i)=a(i-1)/u(i-1);%先求l(i)u(i)=b(i)-c(i-1)*l(i);%再求u(i)%A=LU,Ax=LUx=d,y=Ux,%Ly=d,由于L是下三角矩阵，对角线均为1，所以可求y(i)y(i)=d(i)-l(i)*y(i-1);endx(n)=y(n)/u(n);for i=(n-1):-1:1%Ux=y,由于U是上三角矩阵，所以可求x(i)x(i)=(y(i)-c(i)*x(i+1))/u(i);end2、输入已知参数>>a=[2 2 2 2 2 2 2];>>b=[2 5 5 5 5 5 5 5];>>c=[2 2 2 2 2 2 2];>>d=[220/27 0 0 0 0 0 0 0];3、按定义格式调用函数>>x=zhuiganfa(a,b,c,d)4、输出结果x=[8.147775166909105 -4.073701092835030 2.036477565178471 -1.017492820111148 0.507254485099400 -0.250643392637350 0.119353996493976 -0.047741598597591]第15题【解】1、编写一个程序生成题目条件生成线性方程组A x=b 的系数矩阵A 和右端项量b ，分别定义矩阵A 、B 、a 、b 分别表示系数矩阵，其中1(10.1;,1,2,...,)j ij i i a x x i i j n -==+=或1(,1,2,...,)1ij a i j n i j ==+-分别构成A 、B 对应右端项量分别a 、b 。
数值分析matlab上机作业报告
一、给定向量x ≠0，计算初等反射阵H k 。
1．程序功能：给定向量x ≠0，计算初等反射阵H k 。
2．基本原理： 若()xx x R x ∈=,,, 的分量不全为零，则由12112212122()x (,,,)1()22n T T sign x e x x x x σσσρσσρ-⎧=⎪=+=+⎪⎪⎪==+⎨⎪⎪=-=-⎪⎪⎩u x u uu H I I uuu 确定的镜面反射阵H 使得y e Hx =-=σ；当(1)k n ≤<时，由21/2k ()T 1()()()k 1()()()(())(0,,0,,,,)1()()=()2()nk i i kk nk k k n k T k k Tk k k kk k T k k sign x x x x x x σσρσσσρ=+-⎧=⎪⎪⎪=+∈⎨⎪==+⎪⎪=-⎩∑u R u u u H I u u 有T 121(,,,,,0,,0)n k k k x x x σ-=-∈H x R 算法：(1)输入x ，若x 为零向量，则报错 (2)将x 规范化，{}x x x M ,,,max =如果M ＝0，则报错同时转出停机 否则n i M x x i i ,,2,1, =←(3)计算2x =σ，如果0<1x ，则σσ-= (4))(1x +=σσρ (5)计算1,(1)x σ==+u x u (6)1Tρ-=-H I uu (7)(M ,0,,0)σ=-y(8)按要求输出,结束3．变量说明：x －输入的n维向量；n －n维向量x的维数；M －M是向量x的无穷范数，即x中绝对值最大的一项的绝对值；p －Householder初等变换阵的系数ρ；u －Householder初等变换阵的向量Us －向量x的二范数；x －输入的n维向量；n －n维向量x的维数；p －Householder初等变换阵的系数ρ；u －Householder初等变换阵的向量Uk －数k，H*x=y，使得y的第k+1项到最后项全为零；4．程序代码：(1)function [p,u]=holder2(x)%HOLDER2 给定向量x≠0,计算Householder初等变换阵的p,u%程序功能：函数holder2给定向量x≠0,计算Householder初等变换阵的p,u；%输入：n维向量x；%输出：[p,u]。
数值分析matlab上机报告
不动点迭代：运行结果：代码截图：牛顿迭代：运行结果：代码截图：1）运行结果：代码截图：a二分法运行结果：代码截图：牛顿法：运行结果：代码截图：弦位法：代码截图：第二次上机作业运行结果：代码截图：A欧拉方法运行结果：代码截图：B改进的欧拉方法运行结果：代码截图：C四阶龙格库塔方法：代码截图：第三次上机作业第6章:（7版）P368: 5(e), 8, 9Using1) Gaussian elimination2) Gaussian elimination with partial pivoting.3) Gaussian elimination with scaled partial pivoting and three-digit chopping arithmetic to solvethe following linear systems, and compare the approximations to the actual solution.1.19x1 +2.11x2 − 100x3 + x4 = 1.12,14.2x1 − 0.122x2 + 12.2x3 − x4 = 3.44,100x2 − 99.9x3 + x4 = 2.15,15.3x1 + 0.110x2 − 13.1x3 − x4 = 4.16.Actual solution [0.176, 0.0126,−0.0206,−1.18]解：1）高斯消元法运行结果代码截图：2）部分选主元高斯消元法运行结果：代码截图：3）具有比例因子部分选主元高斯消元法运行结果：代码截图：2.Solve Ax = b using the Crout factorization for tridiagonal systems. Let A be the 10×10 tridiagonal matrix given byaii = 2, ai,i+1 = ai,i−1 = −1, for each i = 2, . . . , 9and a11 = a10,10 = 2, a12 = a10,9 = −1.Let b be the ten-dimensional column vector given byb1 = b10 = 1 and bi = 0, for each i = 2, 3, . . . , 9.运行结果：代码截图：3.1）雅克比迭代运行结果：程序截图：2）GS迭代运行结果代码截图：3）SOR方法运行结果：代码截图：。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
引论习题题目：用二分法求方程x^3-x-1=0在[1,2]内的近似根，要求误差不超过10-3程序：function [root,n]=dichotomy(a,b,eps)% dichotomy：二分法求根函数% f:带求解方程；[a,b]求解区间；n为二分次数if nargin<2 disp('输入变量不够');return;else if nargin>3 disp('输入变量过多');return;else nargin==2eps=1.0e-3;endendif a>b disp('区间输入错误')returnendc=(a+b)/2;if f(c)==0y=cn=1returnendn=0;while abs(b-a)>epsif f(a)*f(c)>0a=c;else b=c;endc=(a+b)/2;n=n+1;endroot=c;end运算结果：定义函数f如下：function [yy] = f(x)%定义了测试函数% x为输入变量yy=x^3-x-1;end在commond window中运行结果如下：>> [root,n]=dichotomy(1,2,1.0e-3)root =1.32471n =10习题一题目一给出概率积分2xxy e dx-=的数据表i 0 1 2 3xi 0.46 0.47 0.48 0.49yi 0.4846555 0.4937452 0.5027498 0.5116683 用二次插值计算，试问：（1）当x=0.472时该积分值等于多少？（2）当x为何值时积分值等于0.5？程序：function [ y,t ] = interpolation( X,Y,x )%拉格朗日插值函数%X为自变量数组，Y为函数值数组，x为插值点，t为插值基函数n=length(X);y=0;t=zeros(1,n);for i=1:nt(i)=1;for j=1:nif i==j continueendt(i)=t(i)*(x-X(j))/(X(i)-X(j));endy=y+t(i)*Y(i);endend运算结果：在commond window中运行结果如下：（1）当x=0.472时该积分值等于多少？>> clear>> X=[0.46 0.47 0.48 0.49];>> Y=[0.4846555 0.4937452 0.5027498 0.5116683];>> x=0.472;>> [ y,t ] = interpolation( X,Y,x )y =0.4956t =-0.0480 0.8640 0.2160 -0.0320（2）当x为何值时积分值等于0.5？>> clear>> X=[0.46 0.47 0.48 0.49];>> Y=[0.4846555 0.4937452 0.5027498 0.5116683];>> y=0.5;>> [ x,t ] = interpolation( Y,X,y )x =0.4769t =-0.0452 0.3356 0.7707 -0.0611题目二构造适合下列数据表的三次插值样条函数S（x）:x -1 0 1 3 y -1 1 3 5 y' 6 1程序：function [ y,m ] =spline( X,Y,x,m1,mn )%样条插值函数%X，Y为输入的自变量、函数值数组；%x为插值点(数组)，y为插值结果(数组)，m为各插值节点的一阶导数值%m1为X(1)的一阶导数;mn为X(n)的一阶导数n=length(X);alpha=zeros(n-2,1);beta=zeros(n-2,1);for i=1:(n-2)alpha(i)=(X(i+1)-X(i))/(X(i+1)-X(i)+X(i+2)-X(i+1));beta(i)=3*((1-alpha(i))*(Y(i+1)-Y(i))/(X(i+1)-X(i))+alpha(i)*(Y(i+2)-Y(i+1))/(X(i+2)-X(i+1) ));endm=zeros(n,1);m(1)=m1;m(n)=mn;A=zeros((n-2),(n-2));B=zeros((n-2),1);for j=1:(n-2)A(j,j)=2;endfor k=2:(n-2)A(k,(k-1))=1-alpha(k);endfor l=1:(n-3)A(l,(l+1))=alpha(l);endB(1)=beta(1)-(1-alpha(1))*m1;B(n-2)=beta(n-2)-alpha(n-2)*mn;for p=2:(n-3)B(p)=beta(p);endm(2:(n-1))=A\B;s=length(x);y=zeros(1,s);for t=1:sr=1;for q=1:nif x(t)>=X(q)&&x(t)<=X(q+1)break;endr=r+1;endxx=(x(t)-X(r))/(X(r+1)-X(r));y(t)=(xx-1)^2*(2*xx+1)*Y(r)+xx^2*(-2*xx+3)*Y(r+1)+(X(r+1)-X(r))*xx*(xx-1)^2*m(r)+( X(r+1)-X(r))*xx^2*(xx-1)*m(r+1);endend运算结果：在commond window中运行结果如下：>> clear>> X=[-1 0 1 3];>> Y=[-1 1 3 5];>> m1=6;mn=1;>> x=-1:0.1:3;>> [ y,m ] =spline( X,Y,x,m1,mn );plot(x,y)求得结果如下（整理之后）：其中红色数据表示插值节点将以上结果绘制成图，如下所示:习题二题目编写一通用型复化辛普生公式，能够对任意长度的等间距离散数据进行积分运算。
要求：命令格式：y=simpson(x,h);y：积分结果x：输入数组h：间距程序：function [ y ] = simpson( x,h )% 一通用型复化辛普森求积公式% y:积分结果；x:输入数组；h:间距if nargin<2 disp('输入参数不够，请准确输入积分数据和步长');return;endn=length(x);% X=zeros(n+1);if n==1 y=x;return;else if mod(n,2)==1% X=zeros(n+1);A=zeros((n+1)/2);B=zeros((n-1)/2);A=x(1:2:(end-2));B=x(2:2:(end-1));y=2*h/6*(x(end)-x(1)+4*sum(B)+2*sum(A));% X(1:n)=x;else mod(n,2)==1X=zeros(n-1);X=x(1:(end-1));m=length(X);A=zeros((m+1)/2);B=zeros((m-1)/2);A=X(1:2:(end-2));B=X(2:2:(end-1));y=2*h/6*(X(end)-X(1)+4*sum(B)+2*sum(A))+(x(n)+x(n-1))/2*h;endendend运算结果：以计算10sin()xI dxx=⎰为例在commond window中运行结果如下：clear>> t=0:(1/8):1;>> x=sin(t)./t;>> h=1/8;>> [ y ] = simpson( x,h )y =0.9461以上积分结果与教材P65所给结果完全一致习题三题目分别用显式和隐式的二阶亚当姆斯方法求解初值问题y'=1-y，y(0)=0，令y(0.2)=0.181，取h=0.2计算y(1.0)。
程序：1、显式亚当姆斯方法：function [ X,Y ] = Adams_shown( x0,y0,y1,h,N )%显式二阶亚当姆斯求解微分方程%x0，y0,y1为输入初值，h为步长，N为计算点数（包括除x0,y0在内的输入初值）%X,Y分别为输出数组,f为y'=f(x,y)X=zeros(N+1,1);Y=zeros(N+1,1);%为X,Y预置空间X=linspace(x0,x0+N*h,N+1);X=X(:);Y(1)=y0;Y(2)=y1;for n=2:Nyy0=f(x0,y0);yy1=f(x0+h,y1);y2=y1+h/2*(3*yy1-yy0);yy2=f(x0+2*h,y2);Y(n+1)=y2;x0=x0+h;y0=y1;y1=y2;endend2、隐式亚当姆斯方法：function [ X,Y ] = Adams_hidden( x0,y0,h,N )%隐式二阶亚当姆斯求解微分方程%x0，y0为输入初值，h为步长，N为计算点数（不包括输入初值）%X,Y分别为输出数组,f为y'=f(x,y)X=zeros(N+1,1);Y=zeros(N+1,1);%为X,Y预置空间X=linspace(x0,x0+N*h,N+1);X=X(:);Y(1)=y0;for n=1:Nyy0=f(x0,y0);% y1=solve('y0+h/2*(f(x0+h,y1)+yy0)-y1','y1');y1=(y0+h/2*(1+yy0))/(1+h/2);yy1=f(x0+h,y1);Y(n+1)=y1;x0=x0+h;y0=y1;endend3、二阶亚当姆斯预报—校正方法：function [ X,Y ] = Adams_shown_hidden( x0,y0,y1,h,N )%二阶亚当姆斯预报-校正求解微分方程%x0，y0,y1为输入初值，h为步长，N为计算点数（包括除x0,y0在内的输入初值）%X,Y分别为输出数组,f为y'=f(x,y)X=zeros(N+1,1);Y=zeros(N+1,1);%为X,Y预置空间X=linspace(x0,x0+N*h,N+1);X=X(:);Y(1)=y0;Y(2)=y1;for n=2:Nyy0=f(x0,y0);yy1=f(x0+h,y1);y2=y1+h/2*(3*yy1-yy0);yy2=f(x0+2*h,y2);y2=y1+h/2*(yy2+yy1);yy2=f(x0+2*h,y2);Y(n+1)=y2;x0=x0+h;y0=y1;y1=y2;endend运算结果：定义函数y'=1-y如下：function [yy] = f(x,y)%定义了测试函数% x，y为输入变量yy=1-y;end在commond window中运行结果如下：>> clear>> x0=0;y0=0;y1=0.181;h=0.2;N=5;>> [ X1,Y1 ] = Adams_shown( x0,y0,y1,h,N );>> [ X2,Y2 ] = Adams_hidden( x0,y0,h,N );>> [ X3,Y3 ] = Adams_shown_hidden( x0,y0,y1,h,N );>> plot(X1,Y1,'r',X2,Y2,'b',X3,Y3,'g');>> legend('显式亚当姆斯','隐式亚当姆斯',,'二阶亚当姆斯校正-预报方法') 最后所得结果如下表所示：根据上述数据绘制成图，如下所示：习题四题目 用牛顿法求下列方程的根，要求计算结果小数点后有4位有效数字（1）30310,2x x x --== （2）20320,1x x x e x --+==程序：function root=NewtonRoot(f,x0,eps)%牛顿法求方程的根% x0为迭代初值%f 为方程表达式：f(x)=0%eps 为迭代精度%root 为满足精度要求的近似根if(nargin==2)eps=1.0e-4;endtol=1;fx=sym(f);dfx=diff(sym(f));%求导数while(tol>eps)fx0=subs(fx,findsym(fx),x0);dfx0=subs(dfx,findsym(dfx),x0); %求该点的导数值 root=x0-fx0/dfx0; %迭代的核心公式 tol=abs(root-x0);x0=root;endformat long g;root=eval(root);end运算结果：在commond window 中运行结果如下：（1）30310,2x x x --== >> clear>> root=NewtonRoot('x^3-3*x-1',2,1.0e-4)root =1.87938524483667（2）20320,1x x x e x --+== >> root=NewtonRoot('x^2-3*x-exp(x)+2',1,1.0e-4)root =0.257530285426349。