matlab数值分析例题
中国科学院大学数值分析MATLAB大作业(计算机学院开设)
7.编程实现题%{在微电机设计计算中需要查磁化曲线表，通常给出的表是磁密B每间隔100高斯磁路每厘米长所需安匝数at的值，下面要解决B从4000至11000区间的查表问题。
为节省计算机存储单元，采用每500高斯存入一个at值，在利用差分公式计算。
从差分表中看到三阶差分近似于0，计算时只需两阶差分。
当4000≤B ≤10500时用牛顿前插公式；当10500≤B≤11000时用牛顿后插公式；试在计算机上编程实现求任一在区间［4000，11000］内的函数值。
%}syms y0y0=[1.38,1.48,1.58,1.69,1.81,1.94,2.10,2.28,2.50,2.76,3.06,3.41,3.83,4.33,4.93];b=input('enter the number of b:');if b>10500&&b<=11000t=(b-11000)/500.0;r=4.93+0.6*t+0.1*t*(t+1)/2;disp(r)elseif b<=10500&&b>=400k=(b-4000)/500;m=floor(k);t=k-m;r=y0(m+1)+(y0(m+2)-y0(m+1))*t+t*(t-1)*(y0(m+3)-2*y0(m+2)+y0(m+1))/2;disp(r)elsedisp('the number is out of consideration')endclear15.用正交多项式(格拉姆－施密特)作最小二乘拟合的程序syms alpha;syms beta;syms a;syms p;x=input('enter the value of x(for example:[1,2,3,4,5]):');y=input('enter the value of y(for example:[1,2,3,4,5]):');w=input('enter the value of weight(for example:[0.1,0.2,0.3,0.4]):');n=input('enter the value of n:');m=length(x);alpha(2)=sum(w*x)/sum(w);a(1)=sum(w*y)/sum(w);p2x= x-alpha(2);alpha(3)=sum(w*x*(subs(p2x,x))^2)/sum(x*( subs(p2x,x))^2);beta(2)=sum(w*( subs(p2x,x))^2)/sum(w);a(2)=sum(w*y*( subs(p2x,x)))/sum(w*( subs(p2x,x))^2);for k=3:m+1alpha(k+1)=sum(w*x*((subs(p,{k,x},[k,x]))^2)/sum(w*( subs(p,{k,x},[k,x]))^2);beta(k)=sum(w*( subs(p,{k,x},[k,x]))^2)/sum(w*( subs(p,{k-1,x},[k-1,x]))^2);a(k)=sum(w*y* subs(p,{k,x},[k,x]))/sum(w*( subs(p,{k,x},[k,x]))^2);p=(x-alpha(k+1))*subs(p,{k,x},[k,x])-beta(k)*subs(p,{k-1,x},[k-1,x]);end,for i=1:n+1F=0;F=F+a(i)* subs(p,{i,t},[i,t]);end。
#数值分析Matlab作业
数值分析编程作业2012年12月第二章14.考虑梯形电阻电路的设计，电路如下：电路中的各个电流{i1,i2,…,i8}须满足下列线性方程组：121232343454565676787822/252025202520252025202520250i i V R i i i i i i i i i i i i i i i i i i i i -=-+-=-+-=-+-=-+-=-+-=-+-=-+=这是一个三对角方程组。
设V=220V ，R=27Ω，运用追赶法，求各段电路的电流量。
Matlab 程序如下：function chase () %追赶法求梯形电路中各段的电流量 a=input('请输入下主对角线向量a='); b=input('请输入主对角线向量b='); c=input('请输入上主对角线向量c='); d=input('请输入右端向量d='); n=input('请输入系数矩阵维数n='); u(1)=b(1); for i=2:nl(i)=a(i)/u(i-1); u(i)=b(i)-c(i-1)*l(i); endy(1)=d(1); for i=2:ny(i)=d(i)-l(i)*y(i-1); endx(n)=y(n)/u(n); i=n-1; while i>0x(i)=(y(i)-c(i)*x(i+1))/u(i); i=i-1;end x输入如下: >> chase请输入下主对角线向量a=[0,-2,-2,-2,-2,-2,-2,-2]; 请输入主对角线向量b=[2,5,5,5,5,5,5,5];请输入上主对角线向量c=[-2,-2,-2,-2,-2,-2,-2,0]; 请输入方程组右端向量d=[220/27,0,0,0,0,0,0,0]; 请输入系数矩阵阶数n=8 运行结果如下：x = 8.1478 4.0737 2.0365 1.0175 0.5073 0.2506 0.1194 0.0477第三章14.试分别用(1)Jacobi 迭代法；(2)Gauss-Seidel 迭代法解线性方程组1234510123412191232721735143231211743511512x x x x x ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥---⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=--⎢⎥⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥---⎣⎦⎣⎦⎣⎦ 迭代初始向量(0)(0,0,0,0,0)T x =。
数值分析matlab程序实例
1，秦九韶算法，求出P(x=3)=2+4x+5x^2+2x^3的值clear all;x=3;n=3;a(1)=2;a(2)=4;a(3)=5;a(4)=2v(1)=a(n+1);for k=2:(n+1);v(k)=x*v(k-1)+a(n-k+2);endp=v(n+1)p=,1132，一次线型插值程序：利用100.121.求115的开方。
clear all;x1=100;x2=121;y1=10;y2=11;x=115;l1=(x-x2)/(x1-x2);l2=(x-x1)/(x2-x1);p1=l1*y1+l2*y2p1=10.71433，分段插值程序，已知为S1（x)为（0，0），（1，1），（2，5）（3,8）上的分段一次插值，求S1（1.5）.clear allx=[0123];y=[0158];n=length(x);a=1.5;for i=2:nif(x(i-1)<=a<x(i));endendH1=y(i-1)+(y(i)-y(i-1))/(x(i)-x(i-1))*(a-x(i-1))H1=3.50004）曲线拟合：用一个5次多项式在区间[0，2π]内逼近函数sin(x)。
clear allX=linspace(0,2*pi,50);Y=sin(X);[P,S]=polyfit(X,Y,5)plot(X,Y,'k*',X,polyval(P,X),'k-')P=-0.00560.0874-0.39460.26850.87970.0102S=R:[6x6double]df:44normr:0.03375)求有理分式的导数clear allP=[3,5,0,-8,1,-5];Q=[10,5,0,0,6,0,0,7,-1,0,-100];[p,q]=polyder(P,Q)6）将以下数据按从小到大排序：4.3 5.7 5.2 1.89.4a=[4.35.75.21.89.4];b(1:100)=0;n=1;b(a*10)=1;for k=1:100a(n)=k/10;if b(k)>0a(n)=k/10;n=n+1;endendaa=1.8000 4.3000 5.2000 5.70009.400010.00007)用二分法求方程x 3-x-1=0在[1,2]内的近似根，要求误差不超过10-3。
数值分析(hilbert矩阵)病态线性方程组的求解matlab程序
（Hilbert 矩阵）病态线性方程组的求解理论分析表明，数值求解病态线性方程组很困难。
考虑求解如下的线性方程组的求解Hx = b ，期中H 是Hilbert 矩阵，()ij n n Hh ，11ij h i j ，i ，j = 1,2，…，n 1.估计矩阵的2条件数和阶数的关系2.对不同的n ，取(1,1,,1)nx K ?，分别用Gauss 消去，Jacobi 迭代，Gauss-seidel 迭代，SOR 迭代和共轭梯度法求解，比较结果。
3.结合计算结果，试讨论病态线性方程组的求解。
第1小题：condition.m %第1小题程序t1=20;%阶数n=20x1=1:t1;y1=1:t1;for i=1:t1H=hilb(i);y1(i)=log(cond(H));endplot(x1,y1);xlabel('阶数n');ylabel('2-条件数的对数（log(cond(H)）');title('2-条件数的对数（log(cond(H)）与阶数n 的关系图');t2=200;%阶数n=200x2=1:t2;y2=1:t2;for i=1:t2H=hilb(i);y2(i)=log(cond(H));endplot(x2,y2);xlabel('阶数n');ylabel('2-条件数的对数（log(cond(H)）');title('2-条件数的对数（log(cond(H)）与阶数n 的关系图');画出Hilbert 矩阵2-条件数的对数和阶数的关系n=200时n=20时从图中可以看出，1）在n小于等于13之前，图像近似直线log（cond（H））~1.519n-1.8332）在n大于13之后，图像趋于平缓，并在一定范围内上下波动，同时随着n的增加稍有上升的趋势第2小题：solve.m%m第2小题主程序N=4000;xGauss=zeros(N,1);xJacobi=zeros(N,1);xnJ=zeros(N,1);xGS=zeros(N,1);xnGS=zeros(N,1);xSOR=zeros(N,1);xnSOR=zeros(N,1);xCG=zeros(N,1);xnCG=zeros(N,1);for n=1:N;x=ones(n,1);t=1.1;%初始值偏差x0=t*x;%迭代初始值e=1.0e-8;%给定的误差A=hilb(n);b=A*x;max=100000000000;%可能最大的迭代次数w=0.5;%SOR迭代的松弛因子G=Gauss(A,b);[J,nJ]=Jacobi(A,b,x0,e,max);[GS,nGS]=G_S(A,b,x0,e,max);[S_R,nS_R]=SOR(A,b,x0,e,max,w);[C_G,nC_G]=CG(A,b,x0,e,max);normG=norm(G'-x);xGauss(n)=normG;normJ=norm(J-x);nJ;xJacobi(n)=normJ;xnJ(n)=nJ;normGS=norm(GS-x);nGS;xGS(n)=normGS;xnGS(n)=nGS;normS_R=norm(S_R-x);nS_R;xSOR(n)=normS_R;xnSOR(n)=nS_R;normC_G=norm(C_G-x);nC_G;xCG(n)=normC_G;xnCG(n)=nC_G;endGauss.m%Gauss消去法function x=Gauss(A,b)n=length(b);l=zeros(n,n);x=zeros(1,n);%消去过程for i=1:n-1for j=i+1:nl(j,i)=A(j,i)/A(i,i);for k=i:nA(j,k)=A(j,k)-l(j,i)*A(i,k);endb(j)=b(j)-l(j,i)*b(i);endend%回代过程x(n)=b(n)/A(n,n);for i=n-1:-1:1c=A(i,:).*x;x(i)=(b(i)-sum(c(i+1:n)))/A(i,i);endJacobi.m%Jacobi迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m 可能最大的迭代次数function [x,n]=Jacobi(A,b,x0,e,m)n=length(A);D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B=D\(L+U);f=D\b;x=B*x0+f;n=1;while norm(x-x0)>ex0=x;x=B*x0+f;n=n+1;if n>mdisp('Jacobi迭代次数过多，迭代可能不收敛');break;endendG_S.m%Gauss-Seidel迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数function [x,n]=G_S(A,b,x0,e,m)n=length(A);D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B=(D-L)\U;f=(D-L)\b;x=B*x0+f;n=1;while norm(x-x0)>ex0=x;x=B*x0+f;n=n+1;if n>mdisp('Gauss-Seidel迭代次数过多，迭代可能不收敛');break;endendSOR.m%SOR超松弛迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数,w松弛因子function [x,n]=SOR(A,b,x0,e,m,w)n=length(A);D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B=(D-w*L)\((1-w)*D+w*U);f=(D-w*L)\b*w;x=B*x0+f;n=1;while norm(x-x0)>ex0=x;x=B*x0+f;n=n+1;if n>mdisp('SOR超松弛迭代次数过多，迭代可能不收敛');break;endendCG.m%CG共轭梯度法,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数function [x,n]=CG(A,b,x0,e,m)r=b-A*x0;p=r;alpha=(r'*r)/(p'*(A*p));x=x0+alpha*p;r1=b-A*x;n=1;while norm(r1)>ebelta=(r1'*r1)/(r'*r);p=r1+belta*p;r=r1;x0=x;alpha=(r'*r)/(p'*(A*p));x=x0+alpha*p;r1=b-A*x;n=n+1;if n>mdisp('CG共轭梯度法迭代次数过多，迭代可能不收敛');break;endend。
牛顿插值法matlab程序例题
牛顿插值法是一种常用的数值分析方法，用于构造一个多项式函数，以便在给定的数据点上进行插值。
这个主题在数学和工程领域中有着广泛的应用，特别是在数据拟合和函数逼近方面。
牛顿插值法的核心思想是通过不断地添加新的数据点来构造一个多项式，并利用已知数据点来确定多项式的系数，从而实现对未知数据点的插值预测。
在Matlab中，实现牛顿插值法并不困难，我们可以利用已有的函数和工具来简化计算过程。
下面，我们将通过一个具体的例题来讲解如何使用Matlab编写牛顿插值法的程序，并分析其结果。
我们需要明确牛顿插值法的数学原理。
给定n个互不相同的节点\(x_0, x_1, ... , x_n\)，以及在这些节点上的函数值\(f(x_0), f(x_1), ... , f(x_n)\)，我们希望构造一个n次插值多项式p(x)，满足p(x_i) = f(x_i)，i=0,1,...,n。
牛顿插值多项式的一般形式为：\[p(x) = a_0 + a_1(x - x_0) + a_2(x - x_0)(x - x_1) + ... + a_n(x -x_0)(x - x_1)...(x - x_{n-1})\]其中，\[a_i\]表示插值多项式的系数。
通过牛顿插值法的迭代过程，可以逐步求解出这些系数，进而得到插值多项式的表达式。
接下来，我们将以一个具体的例题来演示如何在Matlab中实现牛顿插值法。
假设我们有如下的数据点和函数值：\(x = [1, 2, 3, 4]\)\(f(x) = [1, 4, 9, 16]\)我们希望利用这些数据点来构造一个插值多项式，并在给定的区间上进行插值计算。
在Matlab中，可以通过interp1函数来进行插值计算，该函数支持多种插值方法，包括牛顿插值法。
下面是一个简单的Matlab程序示例：```matlabx = [1, 2, 3, 4];y = [1, 4, 9, 16];xi = 2.5;yi = interp1(x, y, xi, 'spline');disp(['在x=',num2str(xi),'处的插值结果为：',num2str(yi)]);```在这段代码中，我们首先定义了给定的数据点x和对应的函数值y，然后利用interp1函数对x=2.5处的插值结果进行计算。
数值分析作业
数值分析作业及答案Chap11、写出下列语句的运行结果。
在MA TLAB 上执行它们以验证所得解答。
a=[1 2 3 ;4 5 6 ]’ b=[9;7;5;3;1] c=b(2:4) d=b(4:-1:1) e=sort(b) f=[3,b ’]解：a=635241 b=13579c=357d=9753 e=97531F=[3 9 7 5 3 1] 3、给定一向量：a=[4 -1 2 -8 4 5 -3 -1 6 -7]写一段程序计算a 中正数的和。
运行程序并显示结果。
解：a=[4 -1 2 -8 4 5 -3 -1 6 -7]; s=0;for i=1:length(a) if a(i)>0s=s+a(i); end end s6、编写一个函数M 文件fun_es(x)，计算如下函数：230.5sin x y e x x =-其中参数可以为标量，也可以为向量。
在MA TLAB 里键入如下命令检验此函数:fun_es(3) fun_es([1 2 3])解：function y=fun_es(x) y=0.5*exp(x/3)-x.^2.*sin(x);chap21、设0x >，x 的相对误差为δ，求L nx 的误差。
解：Lnx-Lnx*=dLnx=dx/x=δ2、设x 的相对误差为2%，求2x 的相对误差。
解：dLnf(x)=xf ’(x)/f(x)dLnx=4%5、计算球体积要使相对误差限为1%，问度量半径R 时允许的相对误差限是多少？解：dLnf(x)=xf ’(x)/f(x)dLnx=3dLnx=1% dLnx=0.33%9、正方形的边长大约为100cm ，应怎样测量才能使其面积误差不超过1cm 2? 解：s=x 2s-s*=2x(x-x*)=1x-x*=1/(2x)=1/200=0.5*10-2 即测量边的误差不超过0.005cm 10、设212S gt =，假定g 是准确的，而对t 的测量有±0.1秒的误差，证明当t 增加时S 的绝对误差增加，而相对误差却减少。
matLAB经典例题及答案
一.对以下数据分别作二次,三次多项式拟合,并画出图形.x=1:16;y=[4,6.4,8,8.4,9.28,9.5,9.7,9.86,10,10.2,10.32,10.42,10.5, 10.55,10.58,10.6];答：程序如下（1）x=(1:16);y=erf(x);p=polyfit(x,y,2);f=polyval(p,x);plot(x,y,x,f);结果p=-0.00100.02020.9096（2）y=[4,6.4,8,8.4,9.28,9.5,9.7,9.86,10,10.2,10.32,10.42,10.5, 10.55,10.58,10.6];y=erf(x);p=polyfit(x,y,3)f=polyval(p,x);plot(x,y,x,f)结果P=0.0002-0.00710.06280.8404二.在[0,4pi]画sin(x),cos(x)(在同一个图象中);其中cos(x)图象用红色小圆圈画.并在函数图上标注“y=sin(x)”,“y=cos(x)”,x轴,y轴,标题为“正弦余弦函数图象”.答：程序如下x=[0:720]*pi/180;plot(x,sin(x),x,cos(x),'ro');x=[2.5;7];y=[0;0];s=['y=sin(x)';'y=cos(x)'];text(x,y,s);xlabel('正弦余弦函数图象'),ylabel('正弦余弦函数图象')图形如下三.选择一个单自由度线性振动系统模型，自定质量、弹簧刚度、阻尼、激振力等一组参数，分别编程（m 文件）计算自由和强迫振动时的响应，并画出振动曲线图。
（要求画出该单自由度线性振动系统模型图）其中质量为m=1000kg,弹性刚度k=48020N/m,阻尼c=1960N.s/m,激振力f(t)=0.阻尼比ζ的程序p=1960/(2*sqrt(48020*1000))求得p=0.1414而p为阻尼比ζ强迫振动时的响应程序g =tf([-101],[48020048020*1.9848020]);bode(g)图形g =tf([001],[0001]);bode(g)振动曲线图程序：函数文件function dx =rigid(t,x)dx =zeros(2,1);dx(1)=x(2);dx(2)=(-48020*x(1)-1960*x(2))/1000;命令文件options =odeset('RelTol',1e-4,'AbsTol',[1e-41e-4]);[T,X]=ode45(@rigid,[012],[11],options);plot(T,X(:,1),'-')其图形如下024681012-6-5-4-3-2-11234单自由度线性强迫振动系统模型图其中质量为m=1000kg,弹性刚度k=48020N/m,阻尼c=1960N.s/m，f(t)=cos(3*pi*t)振动曲线图程序：函数文件function dx=rigid(t,x)dx=zeros(2,1);dx(1)=x(2);dx(2)=(-48020*x(1)-1960*x(2))/1000+cos(3*pi*t);命令文件options=odeset('RelTol',1e-4,'AbsTol',[1e-41e-4]);[T,X]=ode45(@rigid,[020],[11],options);plot(T,X(:,1),'-')力等一组参数，建立Simulink仿真模型框图进行仿真分析。
习题六（Matlab数值计算）课后习题
习题六（Matlab数值计算）课后习题1、利⽤MATLAB 提供的randn 函数⽣成符合正态分布的10×5随机矩阵A,进⾏如下操作：(1)A 各列元素的均值和标准⽅差。
(2)A 的最⼤元素和最⼩元素。
(3)求A 每⾏元素的和以及全部元素之和。
(4)分别对A 的每列元素按升序、每⾏元素按降序排序。
A=randn(10,5);disp('各列元素的均值：');mean(A)disp('各列元素的标准⽅差：');std(A)disp('A 的最⼤元素：');max(max(A))disp('A 的最⼩元素：');min(min(A))disp('A 每⾏元素之和：');sum(A,2)disp('全部元素之和：');sum(sum(A))disp('每列元素按升序：');Y=sort(A)disp('每⾏元素按降序：');Y=sort(A,2,'descend')各列元素的均值：ans =-0.1095 0.1282 -0.2646 0.3030 -0.2464各列元素的标准⽅差：ans =0.9264 1.2631 0.8129 0.8842 1.3151A 的最⼤元素：ans =2.5855A 的最⼩元素：ans =-1.9330A 每⾏元素之和：ans =-2.29701.25450.06615.0489-0.69881.1002-2.9310-2.0595-1.68780.3112全部元素之和：ans =-1.8932每列元素按升序：Y =-1.2141 -1.4916 -1.4224 -1.1658 -1.9330-1.1135 -1.0891 -1.4023 -0.8045 -1.7947-0.8637 -1.0616 -0.7648 -0.2437 -1.1480-0.7697 -0.7423 -0.6156 0.1978 -0.6669-0.2256 0.0326 -0.1961 0.2157 -0.4390-0.0068 0.0859 -0.1924 0.2916 -0.08250.0774 0.5525 -0.1774 0.6966 0.10490.3714 1.1006 0.4882 0.8351 0.18731.1174 1.5442 0.7481 1.4193 0.72231.53262.3505 0.8886 1.5877 2.5855每⾏元素按降序：Y =1.4193 -0.6156 -0.8637 -1.0891 -1.14800.7481 0.2916 0.1049 0.0774 0.03260.7223 0.5525 0.1978 -0.1924 -1.21412.5855 1.5877 1.1006 0.8886 -1.11351.5442 -0.0068 -0.6669 -0.7648 -0.80451.5326 0.6966 0.1873 0.0859 -1.40230.8351 -0.0825 -0.7697 -1.4224 -1.49160.4882 0.3714 -0.2437 -0.7423 -1.93300.2157 -0.1774 -0.2256 -0.4390 -1.06162.3505 1.1174 -0.1961 -1.1658 -1.79472、按要求对指定函数进⾏插值和拟合。
matlab习题及答案
matlab习题及答案《Matlab习题及答案：提升编程技能，解决实际问题》Matlab是一种强大的数值计算和数据分析工具，它被广泛应用于工程、科学和其他领域的计算和模拟。
为了帮助大家提升编程技能，解决实际问题，我们为大家准备了一些Matlab习题及答案，希望能够帮助大家更好地掌握Matlab的使用。
1. 习题一：编写一个Matlab程序，计算斐波那契数列的前20个数字，并将结果打印出来。
答案：下面是一个简单的Matlab程序，用于计算斐波那契数列的前20个数字。
```matlaba = 0;b = 1;fib = zeros(1, 20);fib(1) = a;fib(2) = b;for i = 3:20fib(i) = fib(i-1) + fib(i-2);enddisp(fib);```2. 习题二：编写一个Matlab程序，求解一个二次方程ax^2 + bx + c = 0的根。
答案：下面是一个简单的Matlab程序，用于求解二次方程的根。
```matlaba = 1;b = -3;c = 2;delta = b^2 - 4*a*c;if delta > 0x1 = (-b + sqrt(delta))/(2*a);x2 = (-b - sqrt(delta))/(2*a);disp(['The roots are ', num2str(x1), ' and ', num2str(x2)]);elseif delta == 0x = -b/(2*a);disp(['The root is ', num2str(x)]);elsedisp('The equation has no real roots');end```通过以上两个习题及答案的示例，我们可以看到Matlab的强大功能和灵活性。
通过练习这些习题，我们可以更好地掌握Matlab的基本语法和常用函数，从而在实际问题中更快更准确地解决数值计算和数据分析的挑战。
数值分析大作业(牛顿下山法,拉格朗日法,切比雪夫法)及Matlab程序
课程设计课程名称：数值分析设计题目：学号：姓名：完成时间：2014.11.18题目一： 解线性方程组的直接法 设方程组Ax b =，其中250002511125555111x x x x x x A x x x ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦， 矩阵中10.1(0,1,,5)k x k k =+=，b 由相应的矩阵元素计算，使解向量(1,1,,1)T x =。
(1) A 不变，对b 的元素6b 加一个扰动410-，求解方程组；(2) b 不变，对A 的元素22a 和66a 分别加一个扰动610-，求解方程组； (3) 对上述两种扰动方程组的解做误差分析。
一.数学原理：本计算采用直接法中的列主元高斯消元法，高斯列主元消元法原理如下： 1、设有n 元线性方程组如下：1111n n nn a a a a ⎛⎫ ⎪ ⎪ ⎪⎝⎭1nx x ⎛⎫ ⎪ ⎪ ⎪⎝⎭＝1nb b ⎛⎫ ⎪ ⎪ ⎪⎝⎭2、第一步：如果a11!=0, 令l i1= ai1/a11, I= 2,3,……,n用（-li1）乘第一个方程加到第i 个方程上,得同解方程组：a (1)11 a (1)12 . . . a (1)1nx 1 b (1)1 a (1)21 a (1)22 . . . a (1)2n x 2 b (1)2 . . . . . . . = . a (1)n-11 a (1)n-12 . . a (1)n-1n x n-1 b (1)n-1 a (1)n1 a (1)n2 . . . a (1)nn x n b (1)n简记为：A (2) x = b (2) 其中a (2)ij = a (1)ij – l i1 * a (1)1j , I ,j = 2,3,..,nb (2)I = b (1)I – l i1 * b (1)1 , I = 2,3,...,n 第二步：如果a (2)22 != 0,令l i2= a （2）i2/a （2）22, I= 3,……,n依据同样的原理，对矩阵进行化间（省略），依次下去，直到完成！最后，得到上三角方程组：a(1)11 a(1)12. . . a(1)1nx1b(1)10 a(1)22 . . . a(1)2nx2b(1)2. . . . . . . = .0 0 . . a(n-1)n-1n xn-1b(n-1)n-10 0 . . . a(n)nn xnb(n)n简记为：A(n) x = b(n)最后从方程组的最后一个方程进行回代求解为：Xn = b(n) / a(n)nnXi = ( b(k)k- ∑ a(k)kj x j ) / a(k)kk二．解题过程：1.由题中所给条件可求出b。
(完整word版)matlab数值分析例题
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 —10 0 —10 0 0B=inv（D）*（L+U）% B为雅可比迭代矩阵B =0 0.2500 —0。
25000.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中作出离散点和拟合曲线图。
数值分析Matlab作业龙格库塔欧拉方法解二阶微分方程
% K3 and L3
K4=RK_z(i)+h*L3; L4=rightf_sys1(x(i)+h,RK_y(i)+h*K3,RK_z(i)+h*L3);% K4and L4
RK_y(i+1)=RK_y(i)+1/6*h*(K1+2*K2+2*K3+K4);
RK_z(i+1)=RK_z(i)+1/6*h*(L1+2*L2+2*L3+L4);
end
plot(x,Euler_y,'r-',x,RK_y,'b-');
[y,T]=max(RK_y);
fprintf('角度峰值等于%d',y)%角度的峰值也就是π
运行第三个程序：在一幅图中显示欧拉法和RK4法，随着截断误差的积累，欧拉法产生了较大的误差
h=0.01
h=0.0001，仍然是开始较为稳定，逐渐误差变大
总结：RK4是很好的方法，很稳定，而且四阶是很常用的方法，因为到五阶的时候精度并没有相应提升。通过这两种方法计算出角度峰值y=3.141593，周期是1.777510。
y(0)=0
z(0)=0
精度随着h的减小而更高，因为向前欧拉方法的整体截断误差与h同阶，（因为是用了泰勒公式）所以欧拉方法的稳定区域并不大。
2.RK4-四阶龙格库塔方法
使用四级四阶经典显式Rungkutta公式
稳定性很好，RK4法是四阶方法，每步的误差是h5阶，而总积累误差为h4阶。所以比欧拉稳定。
Matlab作业3(数值分析)答案
Matlab作业3（数值分析）机电工程学院（院、系）专业班组学号姓名实验日期教师评定1.计算多项式乘法(x2+2x+2)(x2+5x+4)。
答：2. （1）将(x-6)(x-3)(x-8)展开为系数多项式的形式。
（2）求解在x=8时多项式(x-1)(x-2) (x-3)(x-4)的值。
答：（1）（2）3. y=sin(x)，x从0到2π，∆x=0.02π，求y的最大值、最小值、均值和标准差。
4.设ｘ＝[0.00.30.8 1.1 1.6 2.3]'，y=[0.500.82 1.14 1.25 1.351.40]'，试求二次多项式拟合系数，并据此计算x1=[0.9 1.2]时对应的y1。
解：x=[0.0 0.3 0.8 1.1 1.62.3]'; %输入变量数据xy=[0.50 0.82 1.14 1.25 1.35 1.40]'; %输入变量数据yp=polyfit(x,y,2) %对x,y用二次多项式拟合,得到系数px1=[0.9 1.2]; %输入点x1y1=polyval(p,x1) %估计x1处对应的y1p =-0.2387 0.9191 0.5318y1 =a) 1.29095.实验数据处理：已知某压力传感器的测试数据如下表p为压力值，u为电压值，试用多项式dcpbpappu+++=23)(来拟合其特性函数，求出a,b,c,d，并把拟合曲线和各个测试数据点画在同一幅图上。
解：>> p=[0.0,1.1,2.1,2.8,4.2,5.0,6.1,6.9,8.1,9.0,9.9];u=[10,11,13,14,17,18,22,24,29,34,39];x=polyfit(p,u,3) %得多项式系数t=linspace(0,10,100);y=polyval(x,t); %求多项式得值plot(p,u,'*',t,y,'r') %画拟和曲线x =0.0195 -0.0412 1.4469 9.8267。
MATLAB解决数值分析问题
1.使用444（x）对数据进行插值，并写出误差分析理论。
建立脚本x1=[0.2 0.4 0.6 0.8 1.0];y1=[0.98 0.92 0.81 0.64 0.38];n=length(y1);c=y1(:);for j=2:n %求差商for i=n:-1:jc(i)=(c(i)-c(i-1))/(x1(i)-x1(i-j+1));endendsyms x df d;df(1)=1;d(1)=y1(1);for i=2:n %求牛顿差值多项式df(i)=df(i-1)*(x-x1(i-1));d(i)=c(i)*df(i);enddisp('4次牛顿插值多项式');P4=vpa(collect((sum(d))),5) %P4即为4次牛顿插值多项式,保留小数点后5位数figureezplot(P4,[0.2,1.08]);输出结果为P4 =- 0.52083*x^4 + 0.83333*x^3 - 1.1042*x^2 + 0.19167*x + 0.98插值余项：R4(x)=f(5)( ξ)/ (5!)* (x - 0.6)*(x - 0.4)*(x - 0.8)*(x - 1)*(x-0.2)新建一个M-filesyms x l;x1=[0.2 0.4 0.6 0.8 1.0];y1=[0.98 0.92 0.81 0.64 0.38];n=length(x1);Ls=sym(0);for i=1:nl=sym(y1(i));for k=1:i-1l=l*(x-x1(k))/(x1(i)-x1(k));endfor k=i+1:nl=l*(x-x1(k))/(x1(i)-x1(k));endLs=Ls+l;endLs=simplify(Ls) %为所求插值多项式Ls(x).figureezplot(Ls,[0.2,1.08]);输出结果为Ls =- (25*x^4)/48 + (5*x^3)/6 - (53*x^2)/48 + (23*x)/120 + 49/50插值余项：R4(x)=f(5)( ξ)/ (5!)* (x - 0.6)*(x - 0.4)*(x - 0.8)*(x - 1)*(x-0.2) 2.试求3次、合曲线，用图示数据曲线及相应的三种拟合曲线。
数值分析实验(4)
页脚内容1实验四 数值积分与数值微分专业班级：信计131班 姓名：段雨博 学号：2013014907 一、实验目的1、熟悉matlab 编程。
2、学习数值积分程序设计算法。
3、通过上机进一步领悟用复合梯形、复合辛普森公式，以及用龙贝格求积方法计算积分的原理。
二、实验题目 P1371、用不同数值方法计算积分049xdx =-⎰。
（1）取不同的步长h .分别用复合梯形及复合辛普森求积计算积分，给出误差中关于h 的函数，并与积分精确值比较两个公式的精度，是否存在一个最小的h ，使得精度不能再被改善？（2）用龙贝格求积计算完成问题（1）。
三、实验原理与理论基础1.1复合梯形公式及其复合辛普森求解[]()()()11101()()222n n n k k k k k h h T f x f x f a f x f b --+==⎡⎤=+=++⎢⎥⎣⎦∑∑误差关于h 的函数：()()212n b a R fh f η-''=-页脚内容2复合辛普森公式：()()()()111/201426n n n k k k k h S f a f x f x f b --+==⎡⎤=+++⎢⎥⎣⎦∑∑误差关于h 的函数：()()441802n n b a h R f I S f η-⎛⎫=-=- ⎪⎝⎭1.2龙贝格求积算法：龙贝格求积公式是梯形法的递推化，也称为逐次分半加速法，它是在梯形公式、辛普森公式和柯特斯公式之间的关系的基础上，构造出一种计算积分的方法，同时它有在不断增加计算量的前提下提高误差的精度的特点。
计算过程如下：（1）取0,k h b a ==-，求：()()()[]()00.,.2hT f a f b k a b =+→⎡⎤⎣⎦令k 1记为区间的二分次数 （2）求梯形值02k b a T -⎛⎫⎪⎝⎭即按递推公式12102122n n n k k h T T f x -+=⎛⎫=+ ⎪⎝⎭∑计算0k T .（3）求加速值，按公式()()()111444141m m k k k mm m m m T T T +--=---逐个求出T 表的地k 行其余各元素()()1,2,,k j j T j k -=（4）若()()001k k T T ε--<（预先给定的精度），则终止计算，并取()()0;1k T I k k ≈+→否则令转（2）继续计算。
数值分析在生活中的应用举例及Matlab实现
一、最小二乘法，用MATLAB实现1. 数值实例下面给定的是乌鲁木齐最近1个月早晨7：00左右（新疆时间）的天气预报所得到的温度，按照数据找出任意次曲线拟合方程和它的图像。
下面用MATLAB编程对上述数据进行最小二乘拟合。
下面用MATLAB编程对上述数据进行最小二乘拟合2、程序代码x=[1:1:30];y=[9,10,11,12,13,14,13,12,11,9,10,11,12,13,14,12,11,10,9,8,7,8,9,11,9,7,6,5,3,1];a1=polyfit(x,y,3) %三次多项式拟合%a2= polyfit(x,y,9) %九次多项式拟合%a3= polyfit(x,y,15) %十五次多项式拟合%b1=polyval(a1,x)b2=polyval(a2,x)b3=polyval(a3,x)r1= sum((y-b1).^2) %三次多项式误差平方和%r2= sum((y-b2).^2) %九次次多项式误差平方和%r3= sum((y-b3).^2) %十五次多项式误差平方和%plot(x,y,'*') %用*画出x,y图像%hold onplot(x,b1, 'r') %用红色线画出x,b1图像%hold onplot(x,b2, 'g') %用绿色线画出x,b2图像%hold onplot(x,b3, 'b:o') %用蓝色o线画出x,b3图像%3、数值结果不同次数多项式拟合误差平方和为：r1=67.6659r2=20.1060r3=3.7952r1、r2、r3分别表示三次、九次、十五次多项式误差平方和。
4、拟合曲线如下图二、 线性方程组的求解（ 高斯-塞德尔迭代算法 ）1、实例： 求解线性方程组（见书P233页）⎪⎪⎩⎪⎪⎨⎧=++=-+=+-3612363311420238321321321x x x x x x x x x 记A x=b, 其中⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--=363320,,12361114238321b x A x x x任取初始值()()Tx0000=，进行迭代。
同济大学数值分析matlab编程
同济⼤学数值分析matlab编程MATLAB 编程题库1.下⾯的数据表近似地满⾜函数21cx bax y ++=，请适当变换成为线性最⼩⼆乘问题，编程求最好的系数c b a ,,，并在同⼀个图上画出所有数据和函数图像.625.0718.0801.0823.0802.0687.0606.0356.0995.0628.0544.0008.0213.0362.0586.0931.0ii y x ----解：>> x=[-0.931 -0.586 -0.362 -0.213 0.008 0.544 0.628 0.995]'; >> y=[0.356 0.606 0.687 0.802 0.823 0.801 0.718 0.625]'; >> A= [x ones(8,1) -x.^2.*y]; >> z=A\y;>> a=z(1); b=z(2); c=z(3); >>xh=[-1:0.1:1]';>>yh=(a.*xh+b)./(1+c.*xh.^2); >>plot(x,y,'r+',xh,yh,'b*')2.若在Matlab ⼯作⽬录下已经有如下两个函数⽂件，写⼀个割线法程序，求出这两个函数精度为1010-的近似根，并写出调⽤⽅式：>> edit gexianfa.mfunction [x iter]=gexianfa(f,x0,x1,tol) iter=0;x=x1;while(abs(feval(f,x))>tol) iter=iter+1;x=x1-feval(f,x1).*(x1-x0)./(feval(f,x1)-feval(f,x0)); x0=x1;x1=x; end>> edit f.m function v=f(x) v=x.*log(x)-1;>> edit g.m function z=g(y) z=y.^5+y-1;>> [x1 iter1]=gexianfa('f',1,3,1e-10) x1 =1.7632 iter1 = 6>> [x2 iter2]=gexianfa('g',0,1,1e-10) x2 =0.7549 iter2 = 83.使⽤GS 迭代求解下述线性代数⽅程组：123123123521242103103x x x x x x x x x ì++=--++=í???-+=??解：>> edit gsdiedai.mfunction [x iter]=gsdiedai(A,x0,b,tol) D=diag(diag(A)); L=D-tril(A); U=D-triu(A); iter=0; x=x0;>> A=[5 2 1;-1 4 2;1 -3 10]; >> b=[-12 10 3]'; >>tol=1e-4; >>x0=[0 0 0]';>> [x iter]=gsdiedai(A,x0,b,tol); >>x x =-3.0910 1.2372 0.9802 >>iter iter = 64.⽤四阶Range-kutta ⽅法求解下述常微分⽅程初值问题（取步长h=0.01）,(1)2x dy y e xy dx y ì??=++?í??=??解：>> edit ksf2.mfunction v=ksf2(x,y)v=y+exp(x)+x.*y; >> a=1;b=2;h=0.01; >> n=(b-a)./h; >> x=[1:0.01:2]; >>y(1)=2;>>for i=2:(n+1)k1=h*ksf2(x(i-1),y(i-1));k2=h*ksf2(x(i-1)+0.5*h,y(i-1)+0.5*k1); k3=h*ksf2(x(i-1)+0.5*h,y(i-1)+0.5*k2); k4=h*ksf2(x(i-1)+h,y(i-1)+k3); y(i)=y(i-1)+(k1+2*k2+2*k3+k4)./6; end >>y调⽤函数⽅法>> edit Rangekutta.mfunction [x y]=Rangekutta(f,a,b,h,y0) x=[a:h:b]; n=(b-a)/h; y(1)=y0; for i=2:(n+1)k1=h*(feval(f,x(i-1),y(i-1)));k2=h*(feval(f,x(i-1)+0.5*h,y(i-1)+0.5*k1)); k3=h*(feval(f,x(i-1)+0.5*h,y(i-1)+0.5*k2)); k4=h*(feval(f,x(i-1)+h,y(i-1)+k3)); y(i)=y(i-1)+ (k1+2*k2+2*k3+k4)./6; end>> [x y]=Rangekutta('ksf2',1,2,0.01,2); >>y5.取0.2h =，请编写Matlab 程序，分别⽤欧拉⽅法、改进欧拉⽅法在12x ≤≤上求解初值问题。
同济大学数值分析matlab编程题汇编.doc
同济大学数值分析matlab编程题汇编.MATLAB编程题库1.下面的数据表近似地满足函数,请适当变换成为线性最小二乘问题,编程求最好的系数,并在同一个图上画出所有数据和函数图像。
解:x=[-x=[:文件一文件二函数v=f(x)v=x . * log(x)-1；函数z=g(y)z=y. y-1；解以下内容:编辑gex AFA。
m函数[x ITER]=gex AFA(f，x0，x1，tol)ITER=0；而(标准(x1-编辑gex AFA。
m函数[x ITER )=gex AFA(f，x0，x1，tol)ITER=0；同时（标准(x1:解以下内容:编辑gsdiedai。
m函数[x iter]=gsdiedai(A，x0，b，tol)D=diag(diag(A))；函数[x iter]=gsdiedai(A，x0，b，tol)D=diag(diag(A))；L=D:编辑ksf 2。
m函数v=ksf2(x，y)v=y exp(x)x . * y；a=1；b=2；h=0.01n=(b-a)./h .x=[1:0.01:2]；y(1)-省略部分-0.5000 1.0000 0.5000-1.0000 1.0000 UU=2.0000 3.0000 4.0000 0-0.5000 7.0000 0 0-1.0000 x=林空间(0，1，11)；x ' ans=0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.90001.0000 x=[1 2 3 4]；y=[6 11 18 27]；p=polyfit(x，y，2)p=1.0000 2.0000 3.0000 diag(1(4，1)，1)ans=0 1 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 012 .编程实现求解满足下列条件的区间[-1，2]上的三次样条函数S(x)，并画出此样条函数的图形: Xi-1 0 1 2 f(Xi)-1 0 1 0f(Xi)' 0-1函数splx=[-1 0 1 2]y=[0-1 0 1 0-1]PP=csape(x，y，'完整')[断点coefs、npolys、ncoefs、dim]=unmkpp(PP)xh=-1:0.133602 if-1=xh=0 yh=coefs(1，1)*(xh 1).3系数(1，2)*(xh 1).2系数(1，3)*(xh 1)系数(1，4)否则，如果0=xh=1 yh=系数(2，1)*(xh).3系数(2，2)*(xh).2系数(2，3)*(xh)系数(2，4)否则，如果1=xh=2 yh=系数(3，1)*(xh-1).3系数(3，2)*(xh-1).2系数(3，3)*(xh-1)系数(3，4)否则返回图(xh，yh，' r ')13 .二分法程序如果nargintol x=(a b)/2 fx=feval(f，x)如果符号(外汇)==符号(fa) a=x fa=fx elseif符号(外汇)==符号b=x FB=FX否则返回结束。
数值分析Matlab作业龙格库塔欧拉方法解二阶微分方程
Matlab 应用使用Euler 和Rungkutta 方法解臂状摆的能量方程背景 单摆是常见的物理模型，为了得到摆角θ的关于时间的函数，来描述单摆运动。
由角动量定理我们知道εJ M =化简得到0sin 22=+θθl g dt d 在一般的应用和计算中，只考虑摆角在5度以内的小摆动，因为可以吧简化为，这样比较容易解。
实际上这是一个解二阶常微分方程的问题。
θ在这里的单摆是一种特别的单摆，具有均匀的质量M 分布在长为2的臂状摆上，使用能量法建立方程W T =h mg ∆=2J 21ω化简得到四阶龙格库塔方法使用四级四阶经典显式Rungkutta公式稳定性很好，RK4法是四阶方法，每步的误差是h5阶，而总积累误差为h4阶。
所以比欧拉稳定。
运行第三个程序：在一幅图中显示欧拉法和RK4法，随着截断误差的积累，欧拉法产生了较大的误差h=0.01h=0.0001，仍然是开始较为稳定，逐渐误差变大总结：RK4是很好的方法，很稳定，而且四阶是很常用的方法，因为到五阶的时候精度并没有相应提升。
通过这两种方法计算出角度峰值y=3.141593，周期是1.777510。
三个程序欧拉法clear;clch=0.00001;a=0;b=25;x=a:h:b;y(1)=0;z(1)=0;for i=1:length(x)-1 % 欧拉y(i+1)=y(i)+h*z(i);z(i+1)=z(i)+h*7.35499*cos(y(i));endplot(x,y,'r*');xlabel('时间');ylabel('角度');A=[x,y];%y(find(y==max(y)))%Num=(find(y==max(y)))[y,T]=max(y);fprintf('角度峰值等于%d',y) %角度的峰值也就是πfprintf('\n')fprintf('周期等于%d',T*h)%周期legend('欧拉');龙格库塔方法先定义函数rightf_sys1.mfunction w=rightf_sys1(x,y,z)w=7.35499*cos(y);clear;clc;%set(0,'RecursionLimit',500)h=0.01;a=0;b=25;x=a:h:b;RK_y(1)=0; %初值%RK_z(1)=0;初值for i=1:length(x)-1K1=RK_z(i); L1=rightf_sys1(x(i),RK_y(i),RK_z(i));%K1 and L1K2=RK_z(i)+0.5*h*L1;L2=rightf_sys1(x(i)+0.5*h,RK_y(i)+0.5*h*K1,RK_z(i)+0.5*h*L1);K3=RK_z(i)+0.5*h*L2;L3=rightf_sys1(x(i)+0.5*h,RK_y(i)+0.5*h*K2,RK_z(i)+0.5*h*L2);K4=RK_z(i)+h*L3;% K4L4=rightf_sys1(x(i)+h,RK_y(i)+h*K3,RK_z(i)+h*L3);and L4RK_y(i+1)=RK_y(i)+1/6*h*(K1+2*K2+2*K3+K4);RK_z(i+1)=RK_z(i)+1/6*h*(L1+2*L2+2*L3+L4);endplot(x,RK_y,'b+');xlabel('Variable x');ylabel('Variable y');A=[x,RK_y];[y,T]=max(RK_y);legend('RK4方法');fprintf('角度峰值等于%d',y) %角度的峰值也就是πfprintf('\n')%周期fprintf('周期等于%d',T*h)两个方法在一起对比使用跟上一个相同的函数rightf_sys1.mclear;clc; %清屏h=0.0001;a=0;b=25;x=a:h:b;Euler_y(1)=0;%欧拉的初值Euler_z(1)=0;RK_y(1)=0;%龙格库塔初值RK_z(1)=0;for i=1:length(x)-1%先是欧拉法Euler_y(i+1)=Euler_y(i)+h*Euler_z(i);Euler_z(i+1)=Euler_z(i)+h*7.35499*cos(Euler_y(i));%龙格库塔K1=RK_z(i); L1=rightf_sys1(x(i),RK_y(i),RK_z(i)); % K1 andL1K2=RK_z(i)+0.5*h*L1;L2=rightf_sys1(x(i)+0.5*h,RK_y(i)+0.5*h*K1,RK_z(i)+0.5*h*L1);% K2 and L2K3=RK_z(i)+0.5*h*L2;L3=rightf_sys1(x(i)+0.5*h,RK_y(i)+0.5*h*K2,RK_z(i)+0.5*h*L2);% K3 and L3K4=RK_z(i)+h*L3; L4=rightf_sys1(x(i)+h,RK_y(i)+h*K3,RK_z(i)+h*L3); K4 and L4RK_y(i+1)=RK_y(i)+1/6*h*(K1+2*K2+2*K3+K4);RK_z(i+1)=RK_z(i)+1/6*h*(L1+2*L2+2*L3+L4);endplot(x,Euler_y,'r-',x,RK_y,'b-');[y,T]=max(RK_y);%角度的峰值也就是πfprintf('角度峰值等于%d',y)fprintf('\n')%周期fprintf('周期等于%d',T*h)xlabel('时间');ylabel('角度');legend('欧拉','RK4');。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
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。