数值分析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
#数值分析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计算实验报告
数值分析第四章外推法计算数值微分MATLAB计算实验报告数值分析MATLAB计算实验报告姓名班级学号⼀、实验名称⽤MATLAB编程实现数值微分的外推法计算。
⼆、实验⽬的1．掌握数值微分和定义和外推法的计算过程；2．了解数值微分外推法的计算⽅法并且编写出与其算法对应的MATLAB程序代码；3．体会利⽤MATLAB软件进⾏数值计算。
三、实验内容⽤外推法计算f(x)=x2e?x在x=0.5的导数。
四、算法描述1.命名函数。
2.如果输⼊未知数少于四个，默认精度10^-33.描述T表矩阵坐标4.依次赋值计算 T表第⼀列5.根据数值微分计算公式求出T表矩阵的值6.若达到精度则运算结束，若未达到循环计算7.输出T表，得出的值就是导数值五、实验结果六、实验结果分析此实验通过MATLAB实现外推法数值微分计算，得到相应的数据，⽅便对数据进⾏分析。
从结果可以看出，当步长h=0.025时⽤中点微分公式只有3位有效数字，外推⼀次达到5位有效数字，外推两次达到9位有效数字。
七、附录（程序）function g=waituifa(fname,x,h,e)if nargin<4,e=1e-3;end;i=1;j=1;G(1,1)=(feval(fname,x+h)-feval(fname,x-h))/(2*h);G(i+1,1)=(feval(fname,x+h/2)-feval(fname,x-h/2))/h;G(i+1,j+1)=(4^j*G(i+1,j)-G(i,j))/(4^j-1);while abs(G(i+1,i+1)-G(i+1,i))>ei=i+1;G(i+1,1)=(feval(fname,x+h/2^i)-feval(fname,x-h/2^i))/(2*h/2^i); for j=1:iG(i+1,j+1)=((4^j)*G(i+1,j)-G(i,j))/(4^j-1);endendGg=G(i+1,i+1);。
数值分析实验— MATLAB实现
数值分析实验——MATLAB实现姓名sumnat学号2013326600000班级13级应用数学2班指导老师2016年1月一、插值：拉格朗日插值 (1)1、代码： (1)2、示例： (1)二、函数逼近：最佳平方逼近 (2)1、代码： (2)2、示例： (2)三、数值积分：非反常积分的Romberg算法 (3)1、代码： (3)2、示例： (4)四、数值微分：5点法 (5)1、代码： (5)2、示例： (6)五、常微分方程：四阶龙格库塔及Adams加速法 (6)1、代码：四阶龙格库塔 (6)2、示例： (7)3、代码：Adams加速法 (7)4、示例： (8)六、方程求根：Aitken 迭代 (8)1、代码： (8)2、示例： (9)七、线性方程组直接法：三角分解 (9)1、代码： (9)2、示例： (10)八、线性方程组迭代法：Jacobi法及G-S法 (11)1、代码：Jacobi法 (11)2、示例： (12)3、代码：G-S法 (12)4、示例： (12)九、矩阵的特征值及特征向量：幂法 (13)1、代码： (13)2、示例： (13)一、插值：拉格朗日插值1、代码：function z=LGIP(x,y)%拉格朗日插值n=size(x);n=n(2);%计算点的个数syms a;u=0;%拉格朗日多项式f=1;%插值基函数for i=1:nfor j=1:nif j==if=f;elsef=f*(a-x(j))/(x(i)-x(j));endendu=u+y(i)*f;f=1;endz=expand(u);%展开2、示例：>> x=1:6;y1=x.^5+3*x.^2-6;y2=sin(x)+sqrt(x);>> f1=LGIP(x,y1)f1 =-6+3*a^2+a^5%可知多项式吻合得很好>> f2=vpa(LGIP(x,y2),3)f2 =.962e-1*a^4+1.38*a+.300*a^2+.504-.436*a^3-.616e-2*a^5二、函数逼近：最佳平方逼近1、代码：function z=BestF(u,a,b,n)%最佳平方逼近，用x^i逼近,n为逼近的次数n=n+1;syms xreal;old=findsym(u);u=subs(u,old,x); %将u中变量替换为xf=sym('');H=sym('');d=sym('');for i=1:n %生成函数系f(1,i)=x^(i-1);endfor i=1:n %生成内积Hfor j=1:nH(i,j)=int(f(1,i)*f(1,j),a,b);endendfor i=1:n %生成内积dd(i,1)=int(f(1,i)*u,a,b);enda=H\d;%解法方程Ha=dz=a'*f';2、示例：>> syms x>> f1=sqrt(x);>> f2=x^5+x^2;>> f3=exp(x);>> a=0 ;b=1;>> BestF(f1,a,b,5)ans =12/143+420/143*x-1120/143*x^2+2016/143*x^3-1800/143*x^4+56/13*x^5>> BestF(f2,a,b,5)ans =x^5+x^2>> BestF(f3,a,b,5)ans =-566826+208524*exp(1)+(16733010-6155730*exp(1))*x+(-115830120+42611520*exp(1))* x^2+(306348840-112699440*exp(1))*x^3+(-342469260+125987400*exp(1))*x^4+(136302012-5 0142708*exp(1))*x^5>> vpa(ans,3)ans =.1e4-.1e6*x-.1e7*x^3+.1e7*x^4三、数值积分：非反常积分的Romberg算法1、代码：function z=IntRom(f,a,b) %Romberg 算法e=1e-10;I{1}=linspace(a,b,2);%1等分I{2}=linspace(a,b,3);%2等分L=setdiff(I{2},I{1});%新得插值点h=b-a;T(1,1)=h/2*sum(subs(f,I{1}));T(2,1)=1/2*T(1,1)+h/2*sum(subs(f,L));T(2,2)=4/3*T(2,1)-1/3*T(1,1);k=2;while abs(T(k,k)-T(k-1,k-1))>e %精度要求k=k+1;I{k}=linspace(a,b,2^(k-1)+1);L=setdiff(I{k},I{k-1});%集合差运算，新得插值点h=h/2;T(k,1)=1/2*T(k-1,1)+h/2*sum(subs(f,L));%梯形for i=2:kT(k,i)=(4^(i-1)/(4^(i-1)-1))*T(k,i-1)-(1/(4^(i-1)-1))*T(k-1,i-1);%加速endEndz=T(k,k);2、示例：>> syms x>> f=x^4;>> a=-4;b=4;>> IntRom(f,a,b)T =1.0e+003 *2.04800000000000 0 0 01.02400000000000 0.68266666666667 0 00.57600000000000 0.42666666666667 0.40960000000000 00.45200000000000 0.41066666666667 0.40960000000000 0.40960000000000ans =4.096000000000000e+002>> vpa((int(f,a,b)-ans),3)ans =0.>> f=exp(x);>> a=0;b=1;>> IntRom(f,a,b)T =Columns 1 through 41.85914091422952 0 0 01.75393109246483 1.71886115187659 0 01.72722190455752 1.71831884192175 1.71828268792476 01.72051859216430 1.71828415469990 1.71828184221844 1.718281828794531.71884112857999 1.71828197405189 1.71828182867536 1.718281828460391.71842166031633 1.71828183756177 1.71828182846243 1.71828182845905Columns 5 through 60 00 00 00 01.71828182845908 01.71828182845905 1.71828182845905ans =1.71828182845905>> vpa((int(f,a,b)-ans),3)ans =0.四、数值微分：5点法1、代码：function z=VDiff(f,x0)%5点法求导数值e=1e-15;h=0.01;for i=0:4x(i+1)=x0+i*h;endy=subs(f,x);m(1)=(1/(12*h))*(-25*y(1)+48*y(2)-36*y(3)+16*y(4)-3*y(5));%5点导数公式h=h/2;for i=0:4x(i+1)=x0+i*h;endy=subs(f,x);m(2)=(1/(12*h))*(-25*y(1)+48*y(2)-36*y(3)+16*y(4)-3*y(5));h=h/2;for i=-0:4x(i+1)=x0+i*h;endy=subs(f,x);m(3)=(1/(12*h))*(-25*y(1)+48*y(2)-36*y(3)+16*y(4)-3*y(5));k=3;while abs(m(k)-m(k-1))<abs(m(k-1)-m(k-2)) & abs(m(k)-m(k-1))>e & (h/10)>0%控制收敛条件及精度要求及h非0h=h/2;k=k+1;for i=0:4x(i+1)=x0+i*h;endy=subs(f,x);m(k)=(1/(12*h))*(-25*y(1)+48*y(2)-36*y(3)+16*y(4)-3*y(5));ende=abs(m(k)-m(k-1));z=[m(k);e];2、示例：>> syms x>> f=exp(x);>> x0=2;>> VDiff(f,x0)ans =7.389056098949710.00000000002558五、常微分方程：四阶龙格库塔及Adams加速法1、代码：四阶龙格库塔function z=RGFour(f,y0,a,b)%4阶龙格库塔，f为函数句柄h=0.01;X=a:h:b;Y(1)=y0;n=size(X);n=n(2);for i=1:n-1K1=f([X(i) Y(i)]);K2=f([X(i)+h/2,Y(i)+h/2*K1]);K3=f([X(i)+h/2,Y(i)+h/2*K2]);K4=f([X(i)+h,Y(i)+h*K3]);Y(i+1)=Y(i)+h/6*(K1+2*K2+2*K3+K4);endz=Y;plot(X,Y);2、示例：>> f=@(x)sin(x(1));>> y0=0;a=0;b=2*pi;>> figure(1);>> RGFour(f,y0,a,b);3、代码：Adams加速法function z=myAdams(f,y0,a,b)h=0.01;p(4)=1;c(4)=1;X=a:h:b;Y(1)=y0;n=size(X);n=n(2);for i=1:3K1=f([X(i) Y(i)]);K2=f([X(i)+h/2,Y(i)+h/2*K1]);K3=f([X(i)+h/2,Y(i)+h/2*K2]);K4=f([X(i)+h,Y(i)+h*K3]);Y(i+1)=Y(i)+h/6*(K1+2*K2+2*K3+K4);endfor i=4:n-1p(i+1)=Y(i)+h/24*(55*f([X(i),Y(i)])-59*f([X(i-1),Y(i-1)])+37*f([X(i-2 ),Y(i-2)])-9*f([X(i-3),Y(i-3)]));m(i+1)=p(i+1)+251/720*(c(i)-p(i));m(i+1)=f([X(i+1),m(i+1)]);c(i+1)=Y(i)+h/24*(9*f([X(i+1),m(i+1)])+19*f([X(i),Y(i)])-5*f([X(i-1), Y(i-1)])+f([X(i-2),Y(i-2)]));Y(i+1)=c(i+1)-19/720*(c(i+1)-p(i+1));endz=Y;plot(X,Y);4、示例：>> f=@(x)exp(x(1));>> myAdams(f,0,0,2*pi);六、方程求根：Aitken 迭代1、代码：function z=myAitken(f,x0);%Aitken 迭代求方程的根e=1e-15;xu1=x0;xv1=subs(f,xu1);xv2=subs(f,xv1);if xv2-2*xv1+xu1==0%防止除0；xu2=xv2;elsexu2=xv2-(xv2-xv1)^2/(xv2-2*xv1+xu1);endwhile abs(xu2-xu1)>e%精度控制xu1=xu2;xv1=subs(f,xu1);xv2=subs(f,xv1);if xv2-2*xv1+xu1==0%防止除0；xu2=xv2;elsexu2=xv2-(xv2-xv1)^2/(xv2-2*xv1+xu1);%Aitken加速公式endendz=xu2;2、示例：>> syms x>> f=cos(x/2)+x;>> x0=3;>> myAitken(f,x0)ans =3.14159265358979>> f=x^2-2+x;>> x0=1;>> myAitken(f,x0)ans =1.41421356237309七、线性方程组直接法：三角分解1、代码：function z=myGuess(A,b);%线性方程组三角分解求根n=size(A);if n~=rank(A)z=['矩阵A线性相关，不符合要求'];return;endn=n(2);L=eye(n);for i=1:n-1for j=i+1:nL(j,i)=A(j,i)/A(i,i);A(j,:)=A(j,:)-L(j,i)*A(i,:);endendU=A;for i=1:n %解Ly=b，得ys=0;for j=1:i-1s=s+y(j)*L(i,j);endy(i)=(b(i)-s)/L(i,i);endfor i=n:-1:1 %解Ux=y，得xs=0;for j=i+1:ns=s+x(j)*U(i,j);endx(i)=(y(i)-s)/U(i,i);endLUz=x';2、示例：>> A=[1 2 3;2 1 5;11 17 21];>> b=[1 3 5]';>> myGuess(A,b)L =1.00000000000000 0 02.00000000000000 1.00000000000000 011.00000000000000 1.66666666666667 1.00000000000000U =1.000000000000002.000000000000003.000000000000000 -3.00000000000000 -1.000000000000000 0 -10.33333333333333ans =-0.06451612903226-0.580645161290320.74193548387097>> t=A\bt =-0.06451612903226-0.580645161290320.74193548387097八、线性方程组迭代法：Jacobi法及G-S法1、代码：Jacobi法function z=myJacobi(A,b)n=size(A);n=n(2);x{1}=zeros(n,1);%初始值e=1e-10;D=diag(diag(A));L=D-tril(A);U=D-triu(A);B=D\(L+U);f=D\b;Q=B'*B;[w,d]=eig(Q);p=max(abs(diag(d))) ; %谱半径if p>=1z=['迭代发散'];return;endx{2}=B*x{1}+f;k=2;while norm(x{k}-x{k-1})>ek=k+1;x{k}=B*x{k-1}+f;endz=x{k};2、示例：>> A=[8 -3 2;4 11 -1;6 3 12];>> b=[20 33 36]';>> myJacobi(A,b)ans =3.000000000013402.000000000012610.999999999992373、代码：G-S法function z=myGS(A,b)n=size(A);n=n(2);x{1}=zeros(n,1);e=1e-10;D=diag(diag(A));L=D-tril(A);U=D-triu(A);B=(D-L)\U;f=(D-L)\b;Q=B'*B;[w,d]=eig(Q);p=max(abs(diag(d))) ; %谱半径if p>=1z=['迭代发散'];return;endx{2}=B*x{1}+f;k=2;while norm(x{k}-x{k-1})>ek=k+1;x{k}=B*x{k-1}+f;endz=x{k};4、示例：>> A=[8 -3 2;4 11 -1;6 3 12];>> b=[20 33 36]';>> myGS(A,b)ans =3.000000000001351.999999999999160.99999999999954九、矩阵的特征值及特征向量：幂法1、代码：function z=myChar(A);%幂法求主特征值及对应的特征向量e=1e-10;n=size(A);n=n(2);v1=ones(n,1);u1=v1;v2=A*u1;a=min(v2);b=max(v2);if abs(a)>abs(b)c=a;elsec=b;endu2=v2/c;%规范化while norm(u2-u1)>eu1=u2;v2=A*u1;a=min(v2);b=max(v2);if abs(a)>abs(b)c=a;elsec=b;endu2=v2/c;%规范化endz{1}=c;z{2}=v2;2、示例：>> A=[8 -3 2;4 11 -1;6 3 12];>> u=myChar(A)u =[14.00000000046956] [3x1 double]>> u{1}ans =14.00000000046956 >> u{2}ans =4.20000000191478 0.93333333198946 14.00000000046956。
牛顿插值法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处的插值结果进行计算。
Matlab与数值分析实验
>> sum
法 2: >>a=1:100; >>sum=a*a’
法 3: >>n=1;
>>sum=0;
>> while n<=100
sum=sum+n*n;
n=n+1;
end
-6-
>> sum * 选择语句
if expression() statements;
[else statements;]
end 例:编写函数文件 demo3 实现 sgn 函数功能
例 2:编写命令文件 demo1 完成以下操作
-4-
建立数组 a=[1,2,3,...,20],b=[1,3,5,...,39],并求 a,b 内积 操作 1) 主窗口点击新建按钮 2) 在弹出的文本编辑窗口添加 a=1:20 b=1:2:39 sum=a*b'
3) 单击保存按钮 将文件命名为 demo1 保存在例 1 新建文件夹中
4) 在 Command Window 中输入 demo1 并回车 例 3:编写函数文件 demo2,返回输入变量的内积
操作:1) 新建 M 文件,编辑如下: function sum=demo2(a,b) sum=a*b';
2) 保存文件在查询目录下,注意不要修改默认名 3) 在 Command Window 中输入
操作:1)新建 M 文件,并编辑如下 function val=demo3(x) if x>0 val=1; elseif x<0 val=-1; else val=0; end 2) 将文件保存在查询目录内 3) >>demo3(0) >>demo3(90) >>demo3(-12)
数值分析大作业(牛顿下山法,拉格朗日法,切比雪夫法)及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。
(完整版)MATLAB数值分析实例
2016-2017 第一学期数值分析上机实验报告姓名： xxx学号： 20162……．学院：土木工程学院导师：………..联系电话：…………..指导老师：………..目录第一题 (1)1.1题目要求 (1)1.2程序编写 (1)1.3计算结果及分析 (2)第二题 (4)2.1题目要求 (4)2.2程序编写 (4)2.3计算结果及分析 (6)第三题 (7)3.1题目要求 (7)3.2程序编写 (7)3.3计算结果及分析 (8)第四题 (9)4.1题目要求 (9)4.2程序编写 (9)4.3计算结果及分析 (10)第五题 (11)5.1题目要求 (11)5.2程序编写 (12)5.3计算结果及分析 (13)第六题 (17)6.1题目要求 (17)6.2程序编写 (17)6.3计算结果及分析 (18)6.4程序改进 (18)第一题选做的是第（1）小问。
1.1题目要求编出不动点迭代法求根的程序；把??3+4??2-10=0写成至少四种??=g(??)的形式，取初值??0=1.5,进行不动点迭代求根，并比较收敛性及收敛速度。
1.2程序编写1.3计算结果及分析①第一种迭代公式：??=??3+4??2+??-10；matlab计算结果如下：（以下为命令行窗口的内容）2；matlab计算结果如下：②第二种迭代公式：??=√(10-??3)/4（以下为命令窗口内容）2；matlab计算结果如下：③第三种迭代公式：??=√10(??+4)?（以下为命令窗口内容）④第四种迭代公式：??=10(??2+4??)?；matlab计算结果如下：（以下为命令窗口内容）上述4种迭代公式，1、4两种由于在x真实值附近|g`(x)|>1,不满足迭代局部收敛条件，所以迭代序列不收敛。
对于2、3两种式子，由于在x真实值附近|g`(x)|<=L<1,满足迭代局部收敛条件,所以迭代序列收敛。
对于2、3两迭代公式，由于L3<L2,所以第3个迭代公式比第2个迭代公式收敛更快。
数值分析中求解非线性方程的MATLAB求解程序(6种)
数值分析中求解非线性方程的MATLAB求解程序（6种）1.求解不动点function [k,p,err,P]=fixpt(g,p0,tol,max1)%求解方程x=g(x) 的近似值，初始值为p0%迭代式为Pn+1=g(Pn)%迭代条件为：在迭代范围内满足|k|<1(根及附近且包含初值)k为斜率P(1)=p0;for k=2:max1P(k)=feval(g,P(k-1));err=abs(P(k)-P(k-1));relerr=err/(abs(P(k))+eps);p=P(k);if (err<tol)|(relerr<tol)break;endendif k==max1disp('超过了最长的迭代次数')endP=P';2.二分法function [c,err,yc]=bisect(f,a,b,delta)%二分法求解非线性方程ya=feval(f,a);yb=feval(f,b);if ya*yb>0break;endmax1=1+round((log(b-a)-log(delta))/log(2));for k=1:max1c=(a+b)/2;yc=feval(f,c);if yc==0a=c;b=c;elseif yb*yc>0b=c;yb=yc;elsea=c;ya=yc;endif b-a<deltabreak;endendc=(a+b)/2;err=abs(b-a);yc=feval(f,c);3.试值法function [c,err,yc]=regula(f,a,b,delta,epsilon,max1) %试值法求解非线性方程%f(a)和飞(b)异号ya=feval(f,a);yb=feval(f,b);if ya*yb>0disp('Note:f(a)*f(b)>0');endfor k=1:max1dx=yb*(b-a)/(yb-ya);c=b-dx;ac=c-a;yc=feval(f,c);if yc==0break;elseif yb*yc>0b=c;yb=yc;elsea=c;ya=yc;enddx=min(abs(dx),ac);if abs(dx)<delta|abs(yc)<epsilonbreak;endendc;err=abs(b-a)/2;yc=feval(f,c);4.求解非线性方程根的近似位置function R=approot(X,epsilon)%求解根近似位置%为了粗估算方程f(x)=0在区间[a,b]的根的位置，%使用等间隔采样点(xk,f(xk))和如下的评定准则：%f(xk-1)与f(xk)符号相反，%或者|f(xk)|足够小且曲线y=f(x)的斜率在%(xk,f(xk))附近改变符号。
数值分析matlab实验报告
数值分析matlab实验报告数值分析 Matlab 实验报告一、实验目的数值分析是研究各种数学问题数值解法的学科，Matlab 则是一款功能强大的科学计算软件。
本次实验旨在通过使用 Matlab 解决一系列数值分析问题，加深对数值分析方法的理解和应用能力，掌握数值计算中的误差分析、数值逼近、数值积分与数值微分等基本概念和方法，并培养运用计算机解决实际数学问题的能力。
二、实验内容（一）误差分析在数值计算中，误差是不可避免的。
通过对给定函数进行计算，分析截断误差和舍入误差的影响。
例如，计算函数＄f(x) ＝＼sin(x)＄在＄x ＝ 05$ 附近的值，比较不同精度下的结果差异。
（二）数值逼近1、多项式插值使用拉格朗日插值法和牛顿插值法对给定的数据点进行插值，得到拟合多项式，并分析其误差。
2、曲线拟合采用最小二乘法对给定的数据进行线性和非线性曲线拟合，如多项式曲线拟合和指数曲线拟合。
（三）数值积分1、牛顿柯特斯公式实现梯形公式、辛普森公式和柯特斯公式，计算给定函数在特定区间上的积分值，并分析误差。
2、高斯求积公式使用高斯勒让德求积公式计算积分，比较其精度与牛顿柯特斯公式的差异。
（四）数值微分利用差商公式计算函数的数值导数，分析步长对结果的影响，探讨如何选择合适的步长以提高精度。
三、实验步骤（一）误差分析1、定义函数｀compute_sin_error` 来计算不同精度下的正弦函数值和误差。
｀｀｀matlabfunction value, error ＝ compute_sin_error(x, precision)true_value ＝ sin(x)；computed_value ＝ vpa(sin(x)， precision)；error ＝ abs(true_value computed_value)；end｀｀｀2、在主程序中调用该函数，分别设置不同的精度进行计算和分析。
（二）数值逼近1、拉格朗日插值法｀｀｀matlabfunction L ＝ lagrange_interpolation(x, y, xi)n ＝ length(x)；L ＝ 0;for i ＝ 1:nli ＝ 1;for j ＝ 1:nif j ～＝ ili ＝ li （xi x(j)）／（x(i) x(j)）；endendL ＝ L ＋ y(i) li;endend｀｀｀2、牛顿插值法｀｀｀matlabfunction N ＝ newton_interpolation(x, y, xi)n ＝ length(x)；％计算差商表D ＝ zeros(n, n)；D(：， 1) ＝ y'；for j ＝ 2:nfor i ＝ j:nD(i, j) ＝（D(i, j 1) D(i 1, j 1)）／（x(i) x(i j ＋ 1)）；endend％计算插值结果N ＝ D(1, 1)；term ＝ 1;for i ＝ 2:nterm ＝ term （xi x(i 1)）；N ＝ N ＋ D(i, i) term;endend｀｀｀3、曲线拟合｀｀｀matlab％线性最小二乘拟合p ＝ polyfit(x, y, 1)；y_fit_linear ＝ polyval(p, x)；％多项式曲线拟合p ＝ polyfit(x, y, n)；％ n 为多项式的次数y_fit_poly ＝ polyval(p, x)；％指数曲线拟合p ＝ fit(x, y, ＇exp1'）；y_fit_exp ＝ p(x)；｀｀｀（三）数值积分1、梯形公式｀｀｀matlabfunction T ＝ trapezoidal_rule(f, a, b, n)h ＝（b a) ／ n;x ＝ a:h:b;y ＝ f(x)；T ＝ h （（y(1) ＋ y(end)）／ 2 ＋ sum(y(2:end 1)））；end｀｀｀2、辛普森公式｀｀｀matlabfunction S ＝ simpson_rule(f, a, b, n)if mod(n, 2) ～＝ 0error(＇n 必须为偶数'）；endh ＝（b a) ／ n;x ＝ a:h:b;y ＝ f(x)；S ＝ h ／ 3 （y(1) ＋ 4 sum(y(2:2:end 1)）＋ 2 sum(y(3:2:end 2)）＋ y(end)）；end｀｀｀3、柯特斯公式｀｀｀matlabfunction C ＝ cotes_rule(f, a, b, n)h ＝（b a) ／ n;x ＝ a:h:b;y ＝ f(x)；w ＝ 7, 32, 12, 32, 7 ／ 90;C ＝ h sum(w y)；end｀｀｀4、高斯勒让德求积公式｀｀｀matlabfunction G ＝ gauss_legendre_integration(f, a, b)x, w ＝ gauss_legendre(5)；％选择适当的节点数t ＝（b a) ／ 2 x ＋（a ＋ b) ／ 2;G ＝（b a) ／ 2 sum(w f(t)）；end｀｀｀（四）数值微分｀｀｀matlabfunction dydx ＝ numerical_derivative(f, x, h)dydx ＝（f(x ＋ h) f(x h)）／（2 h)；end｀｀｀四、实验结果与分析（一）误差分析通过不同精度的计算，发现随着精度的提高，误差逐渐减小，但计算时间也相应增加。
数值分析在生活中的应用举例及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中的实现，同时给出一个具体的例题进行演示。
一、复合形法基本原理复合形法（Nelder-Mead算法）是一种无导数优化方法，它的基本思想是在函数域内随机选择一些点，通过对这些点进行逐步变换、调整，最终找到极小值点。
其算法流程如下：1. 选定n+1个点作为初始点，其中n为被优化函数的自变量个数。
这些点可以随机生成，也可以通过一定的规则生成。
2. 对这些点进行排序，确定最小值点、最大值点和次大值点。
3. 通过一系列变换操作，将最大值点向中间靠拢，同时保持其他点的位置不变。
4. 根据变换后的结果，重新排序，判断是否满足终止条件。
若不满足，继续进行变换操作，直到满足终止条件为止。
终止条件可以是迭代次数、函数值的变化量或某些停止准则的满足等。
二、复合形法MATLAB程序实现复合形法的MATLAB程序可以通过fminsearch函数实现，其语法格式如下：[x,fval,exitflag,output] = fminsearch(fun,x0,options) 其中，fun为被优化的函数，x0为初始点，options为优化选项。
程序返回变量x为最小值点，fval为最小值，exitflag为程序停止状态，output为程序运行信息。
下面是一个简单的例子，用来演示如何使用MATLAB实现复合形法：1. 定义被优化的函数：f(x,y) = (x-2)^2+(y+1)^2+2function z = myfun(x)z = (x(1)-2)^2 + (x(2)+1)^2 + 2;2. 定义初始点和优化选项：x0 = [-1,1];opts = optimset('Display','iter');3. 运行fminsearch函数：[x,fval,exitflag,output] = fminsearch(@myfun,x0,opts);4. 输出结果：x = 1.9999 -0.9999fval = 2.0000exitflag = 1output =iterations: 30funcCount: 58algorithm: 'Nelder-Mead simplex direct search'message: 'Optimization terminated: relative function value decreasing below OPTIONS.TolFun.'最终输出结果表明，通过30次迭代，算法找到了函数的最小值点为x=1.9999，y=-0.9999，函数值为2.0000。
数值分析matlab数值试验
实验一：误差传播及算法稳定性实验1.21、试验程序：function charpt1_2% 误差传播及算法稳定性实验clc;clear all;promps={'请选择递推关系式，若选E1=1/e,En=1-nEn-1,请输入1，若选EN=0，En-1=(1-En)/n,请输入2：'};I=1;while Iresult=inputdlg(promps,'charpt1_2',1,{'1'});Nb=str2num(char(result));if ((Nb~=1)|(Nb~=2))I=0;endend%%%%%%%%%%%%%%%%%I=1;while Iresult=inputdlg('请输入递推步数 n>=1:','charpt1_2',1,{'10'});steps=str2num(char(result));if (steps>0)&(steps==fix(steps)) %% 如果steps大于0且为整数I=0;endend%%%%%%%%%%%%%%%%%result=inputdlg('请输入计算中所采用的有效数字位数n:','charpt1_2',1,{'5'});Sd=str2num(char(result));format long %% 设置显示精度result=zeros(1,steps); %% 存储计算结果err=result; %% 存储计算的绝对误差值func=result; %% 存储用quadl计算的近似值%%%%%%%%%%%%%%%%%%% 用quadl计算积分近似值for n=1:stepsfun=@(x) x.^n.*exp(x-1);func(n)= quadl(fun,0,1);end%%%%%%%%%%%%%%%%%%% 用自定义算法计算if(Nb==1)digits(Sd);result(1)=subs(vpa(1/exp(1)));for n=2:stepsresult(n)=subs(vpa(1-n*result(n-1)));enderr=abs(result-func);elseif(Nb==2)digits(Sd);result(steps)=0;for n=(steps-1):-1:1result(n)=subs(vpa((1-result(n+1))/(n+1)));enderr=abs(result-func);end%%%%%%%%%%%%%%%%%%% 输出结果数值及图像clf;disp('库函数计算值:');disp(sprintf('%e ',func));disp('递推值:');disp(sprintf('%e ',result));disp('误差值:');disp(sprintf('%e ',err));if(Nb==1)plot([1:steps],result,'-rs',[1:steps],func,':k*',[1:steps],err,'-.bo' );elseif(Nb==2)plot([steps:-1:1],result,'-rs',[steps:-1:1],func,':k*',[steps:-1:1],e rr,'-.bo');endxlabel('第n步');ylabel('计算值');legend('自定义算法结果','库函数计算结果','误差值');grid on2、试验结果：选择递推关系式1，递推步数为10，有效数字为5位，计算结果如下：库函数计算值:3.678794e-001 2.642411e-001 2.072766e-001 1.708934e-001 1.455329e-0011.268024e-001 1.123836e-001 1.009323e-001 9.161229e-002 8.387707e-002递推值:3.678800e-001 2.642400e-001 2.072800e-001 1.708800e-001 1.456000e-0011.264000e-001 1.152000e-001 7.840000e-0022.944000e-001 -1.944000e+000误差值:5.588280e-007 1.117662e-006 3.352927e-006 1.341222e-0056.705713e-005 4.023702e-004 2.816427e-003 2.253226e-002 2.027877e-001 2.027877e+00012345678910第n 步计算值选择递推关系式2，递推步数为10，有效数字为5位，计算结果如下： 库函数计算值:3.678794e-001 2.642411e-001 2.072766e-001 1.708934e-001 1.455329e-001 1.268024e-001 1.123836e-001 1.009323e-001 9.161229e-002 8.387707e-002 递推值:3.678800e-001 2.642400e-001 2.072800e-001 1.708900e-001 1.455300e-001 1.267900e-001 1.125000e-001 1.000000e-001 1.000000e-001 0.000000e+000 误差值:5.588280e-007 1.117662e-006 3.352927e-006 3.412224e-006 2.942873e-006 1.237016e-005 1.164270e-004 9.322618e-004 8.387707e-003 8.387707e-002第n 步计算值选择递推关系式1，递推步数为10，有效数字为6位，计算结果如下： 库函数计算值:3.678794e-001 2.642411e-001 2.072766e-001 1.708934e-001 1.455329e-001 1.268024e-001 1.123836e-001 1.009323e-001 9.161229e-002 8.387707e-002 递推值:3.678790e-001 2.642420e-001 2.072740e-001 1.709040e-001 1.454800e-001 1.271200e-001 1.101600e-001 1.187200e-001 -6.848000e-002 1.684800e+000 误差值:4.411720e-007 8.823378e-007 2.647073e-006 1.058778e-0055.294287e-005 3.176298e-004 2.223573e-003 1.778774e-002 1.600923e-001 1.600923e+00012345678910第n 步计算值选择递推关系式2，递推步数为10，有效数字为6位，计算结果如下： 库函数计算值:3.678794e-001 2.642411e-001 2.072766e-001 1.708934e-001 1.455329e-001 1.268024e-001 1.123836e-001 1.009323e-001 9.161229e-002 8.387707e-002 递推值:3.678800e-001 2.642410e-001 2.072770e-001 1.708930e-001 1.455360e-001 1.267860e-001 1.125000e-001 1.000000e-001 1.000000e-001 0.000000e+000 误差值:5.588280e-007 1.176622e-007 3.529274e-007 4.122239e-007 3.057127e-006 1.637016e-005 1.164270e-004 9.322618e-004 8.387707e-003 8.387707e-002第n 步计算值3、结果分析：很明显第二种递推式结果要比第一种好，式1在第七步后有明显误差，而式2在第三步后基本与近似解一致。
数值分析课程设计报告(MATLAB版)
（2）取右端向量 b 的三位有效数字得 b [1.83 1.08 0.783]T ，求方程组的准确 解 X ，并与 X 的数据 [1 1 1]T 作比较 。说明矩阵的病态性。
算法及相应结果： （1）在 MATLAB 命令窗口里输入如下命令： >> H=[1 1/2 1/3;1/2 1/3 1/4;1/3 1/4 1/5]; b=[11/6 13/12 47/60]'; >> x=H\b 回车得到结果为： x = 1.0000 1.0000 1.0000 (2)紧接着在上题基础上继续输入如下命令： >> c=[1.83 1.08 0.783]'; x1=H\c 回车得到如下结果： x1 = 1.0800 0.5400 1.4400
问题分析：考虑由直线段（2 个点）产生第一个图形（5 个点）的过程，设 P 1 和 P5 分别为原始直线段的两个端点。现在需要在直线段的中间依次插入三个点 。显然， P2 位于 P P2 , P3 , P4 产生第一次迭代的图形（图 1-4） 1 点右端直线段的三分 之一处， P4 点绕 P2 旋转 60 度（逆时针方向）而得到的，故可以处理为向量 P2 P4 经正交变换而得到向量 P2 P3 ，形成算法如下： （1） P2 P 1 (P 5 P 1) / 3 ； （2） P4 P 1 2( P 5 P 1) / 3 ； （3） P3 P2 ( P4 P2 ) AT ； 在算法的第三步中，A 为正交矩阵。
运行结果： 0.0884 0.0580 0.0431 0.0343 0.0285 0.0243 0.0212 0.0188 0.0169 0.0154 0.0141 0.0130 0.0120 0.0112 0.0105 0.0099 0.0094 0.0087 0.0092 0.0042 （2）从 I 30 较粗略的估计值出发，我们不妨取 0.01. 源程序：
同济大学数值分析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 ≤≤上求解初值问题。
(完整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编程题汇编.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方法插值法
其中,
n
【注】
x [a, b] w(x) (x x j ) j 0
（1）误差估计
Rn (x)
M n1 (n 1)!
w( x)
M n1
max
x( a ,b )
f
(n1) (x)
（2）余项与 x、M n1 节点的位置、个数 n 有关
（3）当 f (x)是 n 的多项式时Ln (x) f (x) n
M2 2!
(x
x0 )(x
x1 )
其中，
M2
max
x( x0 , x1 )
f
(x)
x
[
6
,
4
]
，所以
R1
(
5
24
)
sin
2!
4 (5
24
)( 5
6 24
)
4
0.0061
2）
抛物插值误差估计.因为
R2 (x)
M3 3!
(x
x0 )(x
x1)(x
x2 )
其中，
M3
max
x( x0 ,x2 )
f (x)
yiynewtonbackwardxyxicos035yi09394数值分析插值法55埃尔米特插值2n12n2数值分析插值法551埃尔米特插值多项式的存在唯一性2n1数值分析插值法数值分析插值法552埃尔米特插值余项553三次埃尔米特插值多项式maxsinxsin1数值分析插值法56561高次插值的病态性质0908原函数150706050405分段线性插值0302543214321055数值分析插值法562分段低次插值方法563分段低次插值余项090807060504030201090807060504030201543214321数值分析插值法57三次样条插值571三次样条插值572三弯矩法数值分析插值法573三次样条插值的误差估计与收敛性58插值运算的matlab函数581一维插值函数interp1yiinterp1xyximethod?linear?yiinterp1xyxilinear?1200135019
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
数值分析（matlab程序）曹德欣曹璎珞第一章绪论数值稳定性程序，计算P11 试验题一积分function try_stableglobal nglobal aa=input('a=');N = 20;I0 = log(1+a)-log(a);I = zeros(N,1);I(1) = -a*I0+1;for k = 2:NI(k) = - a*I(k-1)+1/k;endII = zeros(N,1);if a>=N/(N+1)II(N) = (1+2*a)/(2*a*(a+1)*(N+1));elseII(N) =(1/((a+1)*(N+1))+1/N)/2;endfor k = N:-1:2II(k-1) = ( - II(k) +1/k) / a;endIII = zeros(N,1);for k = 1:Nn = k;III(k) = quadl(@f,0,1);endclcfprintf('\n 算法1结果算法2结果精确值') for k = 1:N,fprintf('\nI(%2.0f) %17.7f %17.7f %17.7f',k,I(k),II(k),III(k)) endfunction y = f(x)global nglobal ay = x.^n./(a+x);return第二章非线性方程求解下面均以方程y=x^4+2*x^2-x-3为例：1、二分法function y=erfen(a,b,esp)format longif nargin<3 esp=1.0e-4;endif fun(a)*fun(b)<0n=1;c=(a+b)/2;while c>espif fun(a)*fun(c)<0b=c;c=(a+b)/2;elseif fun(c)*fun(b)<0a=c;c=(a+b)/2;else y=c; esp=10000;endn=n+1;endy=c;elseif fun(a)==0y=a;elseif fun(b)==0y=b;else disp('these,nay not be a root in the intercal')endnfunction y=fun(x)y=x^4+2*x^2-x-3;2、牛顿法function y=newton(x0)x1=x0-fun(x0)/dfun(x0);n=1;while (abs(x1-x0)>=1.0e-4) & (n<=100000000)x0=x1;x1=x0-fun(x0)/dfun(x0);n=n+1;endy=x1nfunction y=fun(x)y=x^4+2*x^2-x-3; 3、割线法function y=gexian(x0,x1)x2=x1-fun(x1)*(x1-x0)/(fun(x1)-fun(x0)); %根据初始XO 和X1求X2 n=1;while (abs(x1-x0)>=1.0e-4) & (n<=100000000) %判断两个条件截止 x0=x1; %将x1赋给x0 x1=x2; %将x2赋给x1 x2=x1-fun(x1)*(x1-x0)/(fun(x1)-fun(x0)); %迭代运算 n=n+1; end y=x2 nfunction y=fun(x)y=x^4+2*x^2-x-3;第四章题目：推导外推样条公式：⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛--------1232123211223322~~~~~22~n n n n n n n n d d d d M M M Mδμλμλμλδ，并编写程序与Matlab 的Spline 函数结果进行对比，最后调用追赶法解方程组。
解：设()s x 的二阶导数值''()(0,1,,)i i s x M i n ==，因为()s x 在[]1,i i x x +上是三次多项式，所以''()i s x 在[]1,i i x x +上是一次多项式，且可表示为''11()i ii ii i ix x x x s x M M h h ++--=+ （1） 对式（1）进行两次积分可得：33111()()()()()66k k k k k k k k k k kM Ms x x x x x P x x Q x x h h +++=-+-+-+- （2） 由式（2）根据()k k y s x =和11()k k y s x ++=可得系数为：6k k k k k y m h P h =-，116k k k k k y m hQ h ++=- （3） 将式（3）代入式（2）可得()s x 为：3311111()()()()()()()6666k k k k k k k k k k k k k k k k k M M y m h y m hs x x x x x x x x x h h h h +++++=-+-+--+-- 分别对()k k s x 和1()k k s x -进行求导，并根据两者相等整理得：112k k k k k k M M M d μλ-+++= （4）其中，11k k k k h h h μ--=+，1k k k k h h h λ-=+，11116()k k k k k k k k k y y y y d h h h h +-----=-+但是要解给定的方程组，还需要两个另外的条件，而外推样条插值的条件可通过下面推理得出：令式（1）中的x 等于1x ，可得：32112122()M M M h h h h h =+- （5） 1212122()()n n n n n n n n M MM h h h h h -------=+- （6） 令式（4）中2k =，然后将式（5）代入得：1123222(2)(1)h hM M d h h ++-= （7） 令式（4）中1k n =-，然后将式（6）代入得：1112122(2)(1)n n n n n n n h hM M d h h -------++-= （8） 将方程组（4）和式（7）、式（8）联立展开即得题目所求。
按照推导出的外推样条插值公式编程可得如下M 文件（spline_wt.m ） function spline_wt X=[0 1 2 3 4 5];Y=[0 0.5 2 1.5 3.5 1.9]; % 调用自编程序S = spline_w(X,Y); %调用matlab 提供程序 S1=spline(X,Y);fnplt(S,'r',2); % 作图 hold onfnplt(S1,'b',1);hold onplot(X,Y,'or'); % 画上节点title('红线为自编程序曲线，蓝线为自带程序曲线')% pp的第j行表示第j个三次多项式的4系数并写出分段多项式pp = S.coefsP1 = poly2str(pp(1,:),'(x-0)')P2 = poly2str(pp(2,:),'(x-1)')P3 = poly2str(pp(3,:),'(x-2)')ppp=S1.coefsPP1 = poly2str(ppp(1,:),'(x-0)')PP2 = poly2str(ppp(2,:),'(x-1)')PP3 = poly2str(ppp(3,:),'(x-2)')function sp = spline_w(X,Y)n = length(X);h = diff(X);d = diff(Y)./h;d1(2:n-1)=6*diff(d)./(h(1:n-2)+h(2:n-1));mu(2:n-1)=h(1:n-2)./(h(1:n-2)+h(2:n-1));mu(n-1)=1-h(n-2)/h(n-1);la(2:n-1)=1-mu(2:n-1);la(2)=1-h(1)/h(2);% 计算三对角方程组a = mu(3:n-1);b = 2*ones(n-2,1);b(1)=2+h(1)/h(2);b(n-2)=2+h(n-2)/h(n-1);c = la(2:n-2);u(1:n-2) = d1(2:n-1);% 调用追赶法解方程组M(2:n-1) = tridiag(a,b,c,u); M(1)=(1+h(1)/h(2))*M(2)-M(3)*h(1)/h(2); M(n)=(1+h(n-1)/h(n-2))*M(n-1)-M(n-2)*h(n-1)/h(n-2);% 下面计算分段多项式的四个系数S=zeros(n-1,4);for k=0:n-2S(k+1,1)=(M(k+2)-M(k+1))/(6*h(k+1));S(k+1,2)=M(k+1)/2;S(k+1,3)=d(k+1)-h(k+1)*(2*M(k+1)+M(k+2))/6;S(k+1,4)=Y(k+1);endsp = mkpp(X,S);%转换成Matlab 格式% -------- 求解三对角线性方程组的追赶法--------function x=tridiag(a,b,c,d)n = length(d);x = zeros(n,1);for k = 2:nb(k) = b(k) - a(k-1)*c(k-1)/b(k-1);d(k) = d(k) - a(k-1)*d(k-1)/b(k-1);endx(n) = d(n)/b(n);for k = n-1:-1:1x(k) = ( d(k)-c(k)*x(k+1) ) / b(k);end根据所编写的程序可得三次多项式及系数为（包括自编和提供）：pp =-0.9511 3.3533 -1.9022 0-0.9511 0.5000 1.9511 0.50001.7556 -2.3533 0.0978 2.0000-1.5711 2.9133 0.6578 1.5000-1.5711 -1.8000 1.7711 3.5000P1 =-0.95111 (x-0)^3 + 3.3533 (x-0)^2 - 1.9022 (x-0)P2 =-0.95111 (x-1)^3 + 0.5 (x-1)^2 + 1.9511 (x-1) + 0.5P3 =1.7556 (x-2)^3 - 2.3533 (x-2)^2 + 0.097778 (x-2) + 2ppp =-0.9511 3.3533 -1.9022 0-0.9511 0.5000 1.9511 0.50001.7556 -2.3533 0.0978 2.0000-1.5711 2.9133 0.6578 1.5000-1.5711 -1.8000 1.7711 3.5000PP1 =-0.95111 (x-0)^3 + 3.3533 (x-0)^2 - 1.9022 (x-0)PP2 =-0.95111 (x-1)^3 + 0.5 (x-1)^2 + 1.9511 (x-1) + 0.5PP3 =1.7556 (x-2)^3 - 2.3533 (x-2)^2 + 0.097778 (x-2) + 2根据所编写的程序可得生成的图形如图1所示第五章摘要：针对多项式检验的困难性，文中先将多项式线性化为多元线性问题。