matlab数值分析例题

1、在MATLAB中用Jacobi迭代法讨论线性方程组，1231231234748212515xxxxxxxxx（1）给出Jacobi迭代法的迭代方程，并判定Jacobi迭代法求解此方程组是否收敛。（2）若收敛，编程求解该线性方程组。解（1）：A=[4-11;4-81;-215]%线性方程组系数矩阵A=4-114-81-215D=diag(diag(A))D=4000-80005L=-tril(A,-1)%A的下三角矩阵L=000-4002-10U=-triu(A,1)%A的上三角矩阵U=01-100-1000B=inv(D)*(L+U)%B为雅可比迭代矩阵B=00.2500-0.25000.500000.12500.4000-0.20000r=eigs(B,1)%B的谱半径r=0.33471Jacobi迭代法收敛。(2)在matlab上编写程序如下：A=[4-11;4-81;-215];b=[7-2115]';x0=[000]';[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;%迭代次数whilenorm(x-x0)=epsx0=x;x=B*x0+f;k=k+1;if(k=max1)disp('迭代超过300次，方程组可能不收敛');return;endend2、设有某实验数据如下：序号xy序号xy1-3-3.9980.51.37762-2.5-3.3011911.54033-2-2.4161101.51.57074-1.5-1.42931121.58395-1-0.4597122.51.69896-0.50.377581332.01701（1）在MATLAB中作图观察离散点的结构，用多项式拟合的方法拟合一个合适的多项式函数；（2）在MATLAB中作出离散点和拟合曲线图.解（1）：首先观察离散点的结构，matlab中的程序如下，x=[-3-2.5-2-1.5-1-0.500.511.522.53];y=[-3.99-3.3011-2.4161-1.4293-0.45970.3775811.37761.54031.57071.58391.69892.01];plot(x,y,'r*')图形如下：离散点近似如抛物线，所以用二次多项式拟合，所以matlab程序如下：x=[-3-2.5-2-1.5-1-0.500.511.522.53];y=[-3.99-3.3011-2.4161-1.4293-0.45970.3775811.37761.54031.57071.58391.69892.01];s=polyfit(x,y,2);p=poly2str(s,'t')p=-0.22214t^2+1t+0.74384（2）做出离散点与拟合曲线的程序如下：x=[-3-2.5-2-1.5-1-0.500.511.522.53];f=[-3.99-3.3011-2.4161-1.4293-0.45970.3775811.37761.54031.57071.58391.69892.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()Ixxdx（要求积分精度为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.0001endI=0;switchtypecase1,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;case2,n=2;h=(b-a)/2;I1=0;I2=(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b))/h;whileabs(I2-I1)epsn=n+1;h=(b-a)/n;I1=I2;I2=0;fori=0:n-1x=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));endendI=I2;step=n;end4、在MATLAB中用四阶Runge-Kutta法编程求解微分方程初值问题23(03)01dyyxxdxy，并在MATLAB中画图比较方程的解析解与R-K解的结果。解：第1步：首先先把经典RK4子程序编出来如下，用RK4.m保存。RK4.mmatlab程序如下：function[t,y]=RK4(func,t0,tt,y0,N,varagin)%Rk方法计算一阶级微分方程组，%由微分方程知识可以知道，对于高阶微分方程可以化为一阶微分方程组。%初始时刻为t0，结束时为tt，初始时刻函数值为y0%N为步数ifnargin4N=100;end%如果没有输入N的话，那么默认计算步长为（tt-t0）/100y(1,:)=y0(:)';h=(tt-t0)/N;%步长t=t0+[0:N]'*h;fork=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.mmatlab程序如下：%RK4方法的测试函数functionf=RK_fun(t,y)f=-y+t*t+3;第3步：把用RK方法计算的主函数写出，用RK_main.m保存。RK_main.mmatlab程序如下：%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+5yrkyrk=1.00001.52671.96472.38372.83523.35753.97894.72035.59726.62137.8010最后调用主函数得：RK_main得到分析结果如下图：5、在MATLAB中利用牛顿切线迭代法计算非线性方程324100xx在区间[1，2]上的一个根。解，利用matlab写如下程序得：root=newtonqxf('x*x*x+4*x*x-10',1,2,1.0e-6)root=1.3652附牛顿切线法程序如下：functionroot=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*f20)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(dfadfb)root=a-fa/dfa;elseroot=b-fb/dfb;endwhile(toleps)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);endend