数值分析Matlab作业
数值分析matlab作业龙格库塔欧拉方法解二阶微分方程
Matlab 应用使用Euler 和Rungkutta 方法解臂状摆的能量方程背景 单摆是常见的物理模型，为了得到摆角θ的关于时间的函数，来描述单摆运动。
由角动量定理我们知道εJ M =化简得到0sin 22=+θθlgdt d在一般的应用和计算中，只考虑摆角在5度以内的小摆动，因为可以吧简化为θ，这样比较容易解。
实际上这是一个解二阶常微分方程的问题。
在这里的单摆是一种特别的单摆，具有均匀的质量M 分布在长为2的臂状摆上， 使用能量法建立方程WT =hmg ∆=2J 21ω 化简得到θθcos 35499.722=dtd重力加速度取9.806651使用欧拉法令dxdy z =，这样降阶就把二阶常微分方程转化为一阶微分方程组，再利用向前Euler 方法数值求解。
y(i+1)=y(i)+h*z(i);z(i+1)=z(i)+h*7.35499*cos(y(i)); y(0)=0 z(0)=0精度随着h 的减小而更高，因为向前欧拉方法的整体截断误差与h 同阶，（因为是用了泰勒公式）所以欧拉方法的稳定区域并不大。
2.RK4-四阶龙格库塔方法使用四级四阶经典显式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 L1 K2=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;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,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 and L1 K2=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');文- 汉语汉字编辑词条文，wen，从玄从爻。
中国科学院大学数值分析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)
将系数矩阵 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
#数值分析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。
数值分析作业MATLAB
1.用二分法解方程 x-lnx=2 在区间【2 ，4】内的根方法: 二分法算法：f=inline('x-2-log(x)');a=2;b=4;er=b-a; ya=f(a);er0=.00001;while er>er0x0=.5*(a+b);y0=f(x0);if ya*y0<0b=x0;elsea=x0;ya=y0;enddisp([a,b]);er=b-a;k=k+1;end求解结果：>> answer13 43.0000 3.50003.0000 3.25003.1250 3.25003.1250 3.18753.1250 3.15633.1406 3.15633.1406 3.14843.1445 3.1484 3.1445 3.1465 3.1455 3.1465 3.1460 3.1465 3.1460 3.1462 3.1461 3.1462 3.1462 3.14623.1462 3.1462 3.1462 3.1462 3.1462 3.1462 最终结果为: 3.14622.试编写MATLAB 函数实现Newton 插值，要求能输出插值多项式。
对函数141)(2+=x x f 在区间[-5，5]上实现10次多项式插值。
Matlab 程序代码如下：%此函数实现y=1/(1+4*x^2)的n 次Newton 插值，n 由调用函数时指定 %函数输出为插值结果的系数向量（行向量）和插值多项式 算法：function [t y]=func5(n) x0=linspace(-5,5,n+1)'; y0=1./(1.+4.*x0.^2); b=zeros(1,n+1); for i=1:n+1 s=0; for j=1:i t=1; for k=1:iif k~=jt=(x0(j)-x0(k))*t;end;end;s=s+y0(j)/t;end;b(i)=s;end;t=linspace(0,0,n+1);for i=1:ns=linspace(0,0,n+1);s(n+1-i:n+1)=b(i+1).*poly(x0(1:i));t=t+s;end;t(n+1)=t(n+1)+b(1);y=poly2sym(t);10次插值运行结果：[b Y]=func5(10)b =Columns 1 through 4-0.0000 0.0000 0.0027 -0.0000Columns 5 through 8-0.0514 -0.0000 0.3920 -0.0000Columns 9 through 11-1.1433 0.0000 1.0000Y =- (7319042784910035*x^10)/147573952589676412928 + x^9/18446744073709551616 + (256*x^8)/93425 -x^7/1152921504606846976 -(28947735013693*x^6)/562949953421312 -(3*x^5)/72057594037927936 + (36624*x^4)/93425 -(5*x^3)/36028797018963968 -(5148893614132311*x^2)/4503599627370496 +(7*x)/36028797018963968 + 1b为插值多项式系数向量，Y为插值多项式。
数值分析作业-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 。
数值分析作业
数值分析作业及答案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程序
课程设计课程名称：数值分析设计题目：学号：姓名：完成时间：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中作出离散点和拟合曲线图。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
数值分析编程作业2012年12月第二章14.考虑梯形电阻电路的设计，电路如下：电路中的各个电流{i1,i2,…,i8}须满足下列线性方程组：121232343454565676787822/252025202520252025202520250i i V Ri i ii i ii i ii i ii i ii i ii 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;endx输入如下:请输入下主对角线向量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 =。
（1）雅可比迭代法程序如下：function jacobi() %Jacobi 迭代法 a=input('请输入系数矩阵a='); b=input('请输入右端向量b='); x0=input('请输入初始向量x0='); n=input('请输入系数矩阵阶数n='); er=input('请输入允许误差er='); N=input('请输入最大迭代次数N='); for i=1:n for j=1:n if i==jd(i,j)=a(i,j); elsed(i,j)=0; end end endm=eye(5)-d\a; %迭代矩阵 g=d\b;x=m*x0+g; k=1;while k<=N %进行迭代 for i=1:5if max(abs(x(i)-x0(i))) >er x=m*x+g; k=k+1;xreturnendendcontinueendx程序执行如下：>>jacobi请输入系数矩阵a=[10 1 2 3 4;1 9 -1 2 -3;2 -1 7 3 -5;3 2 3 12 -1;4 -3 -5 -1 15] 请输入右端向量b=[12 -27 14 -17 12]'请输入初始向量x0=[0 0 0 0 0]'请输入系数矩阵阶数n=5请输入允许误差er=1.0e-6请输入最大容许迭代次数N=60x =1.0000-2.00003.0000-2.00001.0000（2）高斯－赛德尔迭代法程序如下：function gs_sdl() %gauss-seiddel迭代法a=input('请输入系数矩阵a=');b=input('请输入右端向量b=');x0=input('请输入初始向量x0=');n=input('请输入系数矩阵阶数n=');er=input('请输入允许误差er=');N=input('请输入最大迭代次数N=');for i=1:nfor j=1:nif i<=jl(i,j)=0;elsel(i,j)=-a(i,j);endendendfor i=1:nfor j=1:nif i<ju(i,j)=-a(i,j);elseu(i,j)=0;endendendfor i=1:nfor j=1:nif i==jd(i,j)=a(i,j);elsed(i,j)=0;endendendm=(d-l)\u; %迭代矩阵g=(d-l)\b;x=m*x0+g;k=1;while k<=Nfor i=1:5if max(abs(x(i)-x0(i))) >erx=m*x+g;k=k+1;elsexreturnendendcontinueendx执行结果如下：>> gs_sdl请输入系数矩阵a=[10 1 2 3 4;1 9 -1 2 -3;2 -1 7 3 -5;3 2 3 12 -1;4 -3 -5 -1 15] 请输入右端向量b=[12 -27 14 -17 12]'请输入初始向量x0=[0 0 0 0 0]'请输入系数矩阵阶数n=5请输入允许误差er=1.0e-6请输入最大容许迭代次数N=60x =1.0000-2.00003.0000-2.00001.0000已知如下矩阵，试用幂法求按模最大的特征值与特征向量。
190668430663034236336168147112303628291-⎡⎤⎢⎥-⎢⎥⎢⎥--⎢⎥-⎣⎦Matlab 程序代码如下：function mifa ()A=input('请输入系数矩阵A='); x0=input('请输入初始列向量x0='); n=input('请输入向量维数n='); er=input('请输入允许误差er=');N=input('请输入最大容许迭代次数N='); k=1; mu=0;while k<=N for t=1:nif abs(x0(t))==max(abs(x0)) alfa=x0(t);xb=t; %最大的x0（i ）的下标 end endy=x0./alfa; x0=A*y;lamda=x0(xb); k=k+1; endlamda %按模最大的特征值x0 %按模最大的特征值对应的特征向量 程序执行结果如下： >> mifa请输入系数矩阵A=[190 66 -84 30;66 303 42 -36;336 -168 147 -112;30 -36 28 291] 请输入初始列向量x0=[0 0 0 1]' 请输入向量维数n=4请输入允许误差er=1.0e-6请输入最大容许迭代次数N=100 lamda = 343.0000 x0 =114.3333 343.0000 -0.0000 -171.5002试编写MATLAB函数实现Newton插值，要求能输出插值多项式。
对函数在区间[-5，5]上实现10次多项式插值。
Matlab程序代码如下：%此函数实现y=1/(1+4*x^2)的n次Newton插值，n由调用函数时指定%函数输出为插值结果的系数向量（行向量）和插值多项式function [t y]=func5(n)x0=linspace(-5,5,n+1)';y0=1./(1.+4.*x0.^2);b=zeros(1,n+1);for i=1:n+1s=0;for j=1:it=1;for k=1:iif k~=jt=(x0(j)-x0(k))*t;end;end;s=s+y0(j)/t;end;b(i)=s;end;t=linspace(0,0,n+1);for i=1:ns=linspace(0,0,n+1);s(n+1-i:n+1)=b(i+1).*poly(x0(1:i));t=t+s;end;t(n+1)=t(n+1)+b(1);y=poly2sym(t);10次插值运行结果：[b Y]=func5(10)b =Columns 1 through 4-0.0000 0.0000 0.0027 -0.0000Columns 5 through 8-0.0514 -0.0000 0.3920 -0.0000Columns 9 through 11-1.1433 0.0000 1.0000Y =- (7319042784910035*x^10)/147573952589676412928 + x^9/18446744073709551616 + (256*x^8)/93425 - x^7/1152921504606846976 - (28947735013693*x^6)/562949953421312 - (3*x^5)/72057594037927936 + (36624*x^4)/93425 - (5*x^3)/36028797018963968 - (5148893614132311*x^2)/4503599627370496 + (7*x)/36028797018963968 + 1b为插值多项式系数向量，Y为插值多项式。
插值近似值：x1=linspace(-5,5,101);x=x1(2:100);y=polyval(b,x)y =Columns 1 through 122.70033.99944.3515 4.0974 3.4926 2.7237 1.9211 1.1715 0.5274 0.0154 -0.3571 -0.5960Columns 13 through 24-0.7159 -0.7368 -0.6810 -0.5709 -0.4278 -0.2704 -0.1147 0.0270 0.1458 0.2360 0.2949 0.3227Columns 25 through 360.3217 0.2958 0.2504 0.1915 0.1255 0.0588 -0.0027 -0.0537 -0.0900 -0.1082 -0.1062 -0.0830Columns 37 through 48-0.0390 0.0245 0.1052 0.2000 0.3050 0.4158 0.5280 0.6369 0.7379 0.8269 0.9002 0.9549Columns 49 through 600.9886 1.0000 0.9886 0.9549 0.9002 0.8269 0.7379 0.6369 0.5280 0.4158 0.3050 0.2000Columns 61 through 720.1052 0.0245 -0.0390 -0.0830 -0.1062 -0.1082 -0.0900 -0.0537 -0.0027 0.0588 0.1255 0.1915Columns 73 through 840.2504 0.2958 0.3217 0.3227 0.2949 0.2360 0.1458 0.0270 -0.1147 -0.2704 -0.4278 -0.5709Columns 85 through 96-0.6810 -0.7368 -0.7159 -0.5960 -0.3571 0.0154 0.5274 1.1715 1.9211 2.7237 3.4926 4.0974Columns 97 through 994.3515 3.9994 2.7003绘制原函数和拟合多项式的图形代码：plot(x,1./(1+4.*x.^2))hold allplot(x,y,'r')xlabel('X')ylabel('Y')title('Runge现象')gtext('原函数')gtext('十次牛顿插值多项式')绘制结果：误差计数并绘制误差图：hold offey=1./(1+4.*x.^2)-yey =Columns 1 through 12-2.6900 -3.9887 -4.3403 -4.0857 -3.4804 -2.7109 -1.9077 -1.1575 -0.5128 -0.0000 0.3733 0.6130Columns 13 through 240.7339 0.7558 0.7010 0.5921 0.4502 0.2943 0.1401 0.0000-0.1169 -0.2051 -0.2617 -0.2870Columns 25 through 36-0.2832 -0.2542 -0.2053 -0.1424 -0.0719 -0.0000 0.0674 0.12540.1696 0.1971 0.2062 0.1962Columns 37 through 480.1679 0.1234 0.0660 0.0000 -0.0691 -0.1349 -0.1902 -0.2270-0.2379 -0.2171 -0.1649 -0.0928Columns 49 through 60-0.0271 0 -0.0271 -0.0928 -0.1649 -0.2171 -0.2379 -0.2270-0.1902 -0.1349 -0.0691 0.0000Columns 61 through 720.0660 0.1234 0.1679 0.1962 0.2062 0.1971 0.1696 0.12540.0674 0.0000 -0.0719 -0.1424Columns 73 through 84-0.2053 -0.2542 -0.2832 -0.2870 -0.2617 -0.2051 -0.1169 0.00000.1401 0.2943 0.4502 0.5921 Columns 85 through 960.7010 0.7558 0.7339 0.6130 0.3733 0.0000 -0.5128 -1.1575 -1.9077 -2.7109 -3.4804 -4.0857 Columns 97 through 99-4.3403 -3.9887 -2.6900plot(x,ey) xlabel('X') ylabel('ey')title('Runge 现象误差图')第六章16、钢包问题。