1、 在MATLAB 中用Jacobi 迭代法讨论线性方程组，1231231234748212515x x x x x x x x x -+=⎧⎪-+=-⎨⎪-++=⎩（1）给出Jacobi 迭代法的迭代方程，并判定Jacobi 迭代法求解此方程组是否收敛。

（2）若收敛，编程求解该线性方程组。

解（1）：A=[4 -1 1;4 -8 1;-2 1 5] %线性方程组系数矩阵 A =4 -1 1 4 -8 1 -2 1 5>> D=diag(diag(A)) D =4 0 0 0 -8 0 0 0 5>> L=-tril(A,-1) % A 的下三角矩阵 L =0 0 0 -4 0 0 2 -1 0>> U=-triu(A,1) % A 的上三角矩阵 U =0 1 -1 0 0 -1 0 0 0B=inv(D)*(L+U) % B 为雅可比迭代矩阵 B =0 0.2500 -0.2500 0.5000 0 0.12500.4000 -0.2000 0>> r=eigs(B,1) %B 的谱半径r =0.3347 < 1Jacobi迭代法收敛。

(2) 在matlab上编写程序如下：A=[4 -1 1;4 -8 1;-2 1 5];>> b=[7 -21 15]';>> x0=[0 0 0]';>> [x,k]=jacobi(A,b,x0,1e-7)x =2.00004.00003.0000k =17附jacobi迭代法的matlab程序如下：function [x,k]=jacobi(A,b,x0,eps)% 采用Jacobi迭代法求Ax=b的解% A为系数矩阵% b为常数向量% x0为迭代初始向量% eps为解的精度控制max1= 300; %默认最多迭代300，超过300次给出警告D=diag(diag(A)); %求A的对角矩阵L=-tril(A,-1); %求A的下三角阵U=-triu(A,1); %求A的上三角阵B=D\(L+U);f=D\b;x=B*x0+f;k=1; %迭代次数while norm(x-x0)>=epsx0=x;x=B*x0+f;k=k+1;if(k>=max1)disp('迭代超过300次，方程组可能不收敛');return;endend2、设有某实验数据如下：（1）在MATLAB中作图观察离散点的结构，用多项式拟合的方法拟合一个合适的多项式函数；（2）在MATLAB中作出离散点和拟合曲线图.解（1）：首先观察离散点的结构，matlab中的程序如下，x=[-3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 ];>> y=[-3.99 -3.3011 -2.4161 -1.4293 -0.4597 0.37758 1 1.3776 1.5403 1.5707 1.5839 1.6989 2.01]; >> plot(x,y,'r*')图形如下：离散点近似如抛物线，所以用二次多项式拟合，所以matlab程序如下：x=[-3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 ];>> y=[-3.99 -3.3011 -2.4161 -1.4293 -0.4597 0.37758 1 1.3776 1.5403 1.5707 1.5839 1.6989 2.01]; >> s=polyfit(x,y,2);>> p=poly2str(s,'t')p =-0.22214 t^2 + 1 t + 0.74384（2）做出离散点与拟合曲线的程序如下：x=[-3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 ];>> f=[-3.99 -3.3011 -2.4161 -1.4293 -0.4597 0.37758 1 1.3776 1.5403 1.5707 1.5839 1.6989 2.01]; >> p=polyfit(x,f,2); >> y=polyval(p,x); >> plot(x,f,'r+',x,y,'k'); >> xlabel('x'); >> ylabel('y'); >> title('拟合');>> gtext(datestr(today))得出的离散点与拟合曲线图像如下：3、在MATLAB 中用复合Simpson 公式编程计算 221()I x x dx -=--⎰（要求积分精度为510- ）解：matlab 程序如下：[I,step] = jfSimpson('-x-x^2',-1,2,2,1.0e-5) I =-4.5000step =4附复合Simpson 程序如下：function [I,step] = jfSimpson(f,a,b,type,eps) %type = 1 辛普森公式 %type = 2 复合辛普森公式if (type==2 && nargin==4)eps=1.0e-4; %缺省精度为0.0001 endI=0;switch type case 1,I=((b-a)/6)*(subs(sym(f),findsym(sym(f)),a)+... 4*subs(sym(f),findsym(sym(f)),(a+b)/2)+... subs(sym(f),findsym(sym(f)),b)); step=1;case 2, n=2;h=(b-a)/2; I1=0;I2=(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b))/h; while abs(I2-I1)>eps n=n+1;h=(b-a)/n; I1=I2; I2=0;for i=0:n-1 x=a+h*i; x1=x+h;I2=I2+(h/6)*(subs(sym(f),findsym(sym(f)),x)+... 4*subs(sym(f),findsym(sym(f)),(x+x1)/2)+... subs(sym(f),findsym(sym(f)),x1)); end end I=I2; step=n; end4、在MATLAB 中用四阶Runge-Kutta 法编程求解微分方程初值问题()23(03)01dy y x x dx y ⎧=-++≤≤⎪⎨⎪=⎩， 并在MATLAB 中画图比较方程的解析解与R-K 解的结果。

解：第1步：首先先把经典RK4子程序编出来如下，用RK4.m 保存。

RK4.m matlab 程序如下：function [t,y] = RK4(func,t0,tt,y0,N,varagin)% Rk方法计算一阶级微分方程组，% 由微分方程知识可以知道，对于高阶微分方程可以化为一阶微分方程组。

% 初始时刻为t0，结束时为tt，初始时刻函数值为y0% N为步数if nargin<4N=100;end% 如果没有输入N的话，那么默认计算步长为（tt-t0）/100y(1,:) = y0(:)';h = (tt-t0)/N; %步长t = t0+[0:N]'*h;for k = 1:Nf1 = h*feval(func,t(k),y(k,:));f1 = f1(:)';% RK中的k1f2 = h*feval(func,t(k) + h/2,y(k,:) + f1/2);f2 = f2(:)';% RK中的k2f3 = h*feval(func,t(k) + h/2,y(k,:) + f2/2);f3 = f3(:)';% RK中的k3f4 = h*feval(func,t(k) + h,y(k,:) + f3);f4 = f4(:)';% RK中的k4y(k + 1,:) = y(k,:) + (f1 + 2*(f2 + f3) + f4)/6;% 所求结果End第2步：然后把微分方程的表达式用一个调试程序RK_fun表示，用RK_fun.m 保存。

RK_fun.m matlab程序如下：%RK4方法的测试函数function f=RK_fun(t,y)f=-y+t*t+3;第3步：把用RK方法计算的主函数写出，用RK_main.m 保存。

RK_main.m matlab程序如下：%rk方法的主函数x0=0;xt=3;y0=1;%符号计算给出分析解y=dsolve('Dy=-y+t*t+3','y(0)=1','t');%RK4方法给出计算结果[x,yrk] = RK4('RK_fun',x0,xt,y0,10);%对应于真解的离散数据yreal=subs(y,x);tol=yreal-yrk; %RK4方法与分析解的误差plot(x,yreal,x,yrk,'+',x,tol,'*')legend('分析解','RK4计算结果','两者误差')yrk;[x,yrk,tol];%yrk为所要计算的结果第4步matlab中输入y; yrk得：y =t^2 - 4/exp(t) - 2*t + 5>> yrkyrk =1.00001.52671.96472.38372.83523.3575 3.97894.72035.59726.62137.8010最后调用主函数得：>> RK_main得到分析结果如下图：5、在MATLAB 中利用牛顿切线迭代法计算非线性方程 324100x x +-= 在区间[1，2]上的一个根。

解，利用matlab 写如下程序得：root=newtonqxf('x*x*x+4*x*x-10',1,2,1.0e-6) root =1.3652附牛顿切线法程序如下：function root=newtonqxf(f,a,b,eps)% 牛顿法求函数 f在区间[a,b]上的一个零点% f 为函数名% a 为区间左端点% b 为区间右端点% eps 为根的精度% root 为求出的函数零点if(nargin==3)eps=1.0e-6;endf1=subs(sym(f),findsym(sym(f)),a);f2=subs(sym(f),findsym(sym(f)),b);if(f1==0)root=a;endif(f2==0)root=b;endif(f1*f2>0)disp('两端点函数值乘积大于0!');return;elsetol=1;fun=diff(sym(f));fa=subs(sym(f),findsym(sym(f)),a);fb=subs(sym(f),findsym(sym(f)),b);dfa=subs(sym(fun),findsym(sym(fun)),a);dfb=subs(sym(fun),findsym(sym(fun)),b);if(dfa>dfb)root=a-fa/dfa;elseroot=b-fb/dfb;endwhile(tol>eps)r1=root;fx=subs(sym(f),findsym(sym(f)),r1);dfx=subs(sym(fun),findsym(sym(fun)),r1); root=r1-fx/dfx;tol=abs(root-r1);end end。