数值分析实验作业matlab编程课题八

合集下载

数值分析上机作业(MATLAB)

数值分析上机作业(MATLAB)
代矩阵。根据迭代矩阵的不同算法,可分为雅各比迭代方法和高斯-赛德尔方法。 (a)雅各比算法
将系数矩阵 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实验报告

数值分析matlab实验报告

数值分析matlab实验报告数值分析MATLAB实验报告引言:数值分析是一门研究利用计算机进行数值计算和解决数学问题的学科。

它在科学计算、工程技术、金融等领域中有着广泛的应用。

本实验旨在通过使用MATLAB软件,探索数值分析的基本概念和方法,并通过实际案例来验证其有效性。

一、插值与拟合插值和拟合是数值分析中常用的处理数据的方法。

插值是通过已知数据点之间的函数关系,来估计未知数据点的值。

拟合则是通过一个函数来逼近一组数据点的分布。

在MATLAB中,我们可以使用interp1函数进行插值计算。

例如,给定一组离散的数据点,我们可以使用线性插值、多项式插值或样条插值等方法,来估计在两个数据点之间的未知数据点的值。

拟合则可以使用polyfit函数来实现。

例如,给定一组数据点,我们可以通过最小二乘法拟合出一个多项式函数,来逼近这组数据的分布。

二、数值积分数值积分是数值分析中用于计算函数定积分的方法。

在实际问题中,往往无法通过解析的方式求得一个函数的积分。

这时,我们可以使用数值积分的方法来近似计算。

在MATLAB中,我们可以使用quad函数进行数值积分。

例如,给定一个函数和积分区间,我们可以使用quad函数来计算出该函数在给定区间上的定积分值。

quad函数使用自适应的方法,可以在给定的误差限下,自动调整步长,以保证积分结果的精度。

三、常微分方程数值解常微分方程数值解是数值分析中研究微分方程数值解法的一部分。

在科学和工程中,我们经常遇到各种各样的微分方程问题。

而解析求解微分方程往往是困难的,甚至是不可能的。

因此,我们需要使用数值方法来近似求解微分方程。

在MATLAB中,我们可以使用ode45函数进行常微分方程数值解。

例如,给定一个微分方程和初始条件,我们可以使用ode45函数来计算出在给定时间范围内的解。

ode45函数使用龙格-库塔方法,可以在给定的误差限下,自动调整步长,以保证数值解的精度。

结论:本实验通过使用MATLAB软件,探索了数值分析的基本概念和方法,并通过实际案例验证了其有效性。

数值分析作业(完整版)

数值分析作业(完整版)

的逆阵 A ,用左除命令 A \ E 检验你的结果。
clc clear close all A=[1 1 1 1 1;1 2 3 4 5;1 3 6 10 15;1 4 10 20 35;1 5 15 35 70]; fprintf('对上述矩阵进行列主元素分解:\n') for i=1:1:r-1 [mx,ro]=max(abs(A(i:r,i))); % 寻找a阵第i列的最大值 [A(i,:),A(ro+i-1,:)]=exchange(A(i,:),A(ro+i-1,:)); % 进行行与行交换 for j=i+1:1:r A(j,:)=A(j,:)-A(j,i)/A(i,i)*A(i,:); end A End %--矩阵A的逆阵 A1=inv(A) %--左除验证 E=eye(5); A2=A\E % 5x5单位阵 % A阵的逆矩阵 % 输出每次交换后的A
第一章
1、计算积分 I n
Code: clc clear close all n=9; %--梯形积分法 x=0:0.01:1; y=(x.^n).*exp(x-1); In = trapz(x,y); In2=vpa(In,6) % 6位有效数字 %--高精度积分法 F = @(x1)(x1.^n).*exp(x1-1); s = quad(F,0,1); s1=vpa(s,6)
0
0, 0, 0, 0, 0 。
T
if abs(er(:,i-1))<=e fprintf('在迭代 %d 次之后,满足精度要求,x向量的值如下:\n',i); fprintf('x1=%.5f, x2=%.5f, x3=%.5f, x4=%.5f, x5=%.5f\n',x(1,i),x(2,i),x(3,i),x(4,i),x(5,i)); break end end %--绘图 figure(1) plot(1:1:i,x(1,:),'b',1:1:i,x(2,:),'k',1:1:i,x(3,:),'g',1:1:i,x(4,:), 'r',1:1:i,x(5,:),'c') legend('x1','x2','x3','x4','x5') grid on title('Jacobi迭代法——x值随迭代次数变化曲线') figure(2) plot(1:1:i-1,er(1,:),'b',1:1:i-1,er(2,:),'k',1:1:i-1,er(3,:),'g',1:1: i-1,er(4,:),'r',1:1:i-1,er(5,:),'c') legend('△x1','△x2','△x3','△x4','△x5') grid on title('Jacobi迭代法——△x值随迭代次数变化曲线') %% fprintf('\n-------------Gauss-Seidel迭代法---------------------\n'); U=-(A1-D); L=-(A2-D); DL_1=inv(D-L); M1=DL_1*U; b2=DL_1*b; x1(:,1)=M1*x0+b2; for j=2:1:100 x1(:,j)=M1*x1(:,j-1)+b2; er1(:,j-1)=x1(:,j)-x1(:,j-1); if abs(er1(:,j-1))<=e fprintf('在迭代 %d 次之后,满足精度要求,x向量的值如下:\n',j); fprintf('x1=%.5f, x2=%.5f, x3=%.5f, x4=%.5f, x5=%.5f\n',x1(1,j),x1(2,j),x1(3,j),x1(4,j),x1(5,j)); break end end %--绘图 figure(3) plot(1:1:j,x1(1,:),'b',1:1:j,x1(2,:),'k',1:1:j,x1(3,:),'g',1:1:j,x1(4 ,:),'r',1:1:j,x1(5,:),'c') legend('x1','x2','x3','x4','x5')

数值分析MATLAB实验报告

数值分析MATLAB实验报告

实验 2.1 多项式插值的震荡现象问题提出:考虑在一个固定的区间上用插值逼近一个函数。

显然Lagrange 插值中使用的节点越多,插值多项式的次数越高,我们自然关心插值多项式的次数增加时,)(x L n 是否也更加靠近被逼近的函数。

Runge 给出的一个例子是极著名并富有启发性的。

设区间[-1,1]上函数.2511)(2xx f +=实验内容:考虑空间[-1,1]的一个等距划分,分点为 nix i 21+-=, =i 0,1,2 ...,n , 则拉格朗日插值多项式为)(2511)(02x l xx L ini n ∑=+=. 其中,n i x l i ,...,2,1,0),(=是n 次Lagrange 插值基函数。

实验要求:(1)选择不断增大的分点数目,...,3,2=n 画出原函数)(x f 及插值多项式函数)(x L n 在[-1,1]上的图像,比较并分析实验结果。

(2)选择其他的函数,例如定义在区间[-5,5]上的函数 ,1)(4xxx h +=)arctan()(x x g =, 重复上述的实验看其结果如何。

首先编写拉格朗日插值函数的Matlab 实现: Matlab 程序为:function y=lagrange(x0,y0,x) %Lagrange 插值 n=length(x0); m=length(x); for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; for j=1:n if(j~=k)p=p*(z-x0(j))/(x0(k)-x0(j)); end ends=s+p*y0(k); endy(i)=s; end(1)当函数为.2511)(2xx f +=时, Matlab 程序为:x=linspace(-1,1,100); y=1./(1+25*x.^2); plot(x,y) hold on; for i=2:2:10x0=linspace(-1,1,i+1); y0=1./(1+25*x0.^2); y=laglanri(x0,y0,x); plot(x,y,'r--') hold on end运行结果:结果分析:从图上看到在区间[-1,1]的两端点附近,随着插值点数的增加,插值函数)(x L n 与)(x f 偏离的越远,而且出现了振荡现象。

数值分析作业MATLAB

数值分析作业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实现

数值分析实验——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实验报告《数值分析MATLAB实验报告》摘要:本实验报告基于MATLAB软件进行了数值分析实验,通过对不同数学问题的数值计算和分析,验证了数值分析方法的有效性和准确性。

实验结果表明,MATLAB在数值分析领域具有较高的应用价值和实用性。

一、引言数值分析是一门研究利用计算机进行数值计算和分析的学科,其应用范围涵盖了数学、物理、工程等多个领域。

MATLAB是一种常用的数值计算软件,具有强大的数值分析功能,能够进行高效、准确的数值计算和分析,因此在科学研究和工程实践中得到了广泛的应用。

二、实验目的本实验旨在通过MATLAB软件对数值分析方法进行实验验证,探究其在不同数学问题上的应用效果和准确性,为数值分析方法的实际应用提供参考和指导。

三、实验内容1. 利用MATLAB进行方程求解实验在该实验中,利用MATLAB对给定的方程进行求解,比较数值解和解析解的差异,验证数值解的准确性和可靠性。

2. 利用MATLAB进行数值积分实验通过MATLAB对给定函数进行数值积分,比较数值积分结果和解析积分结果,验证数值积分的精度和稳定性。

3. 利用MATLAB进行常微分方程数值解实验通过MATLAB对给定的常微分方程进行数值解,比较数值解和解析解的差异,验证数值解的准确性和可靠性。

四、实验结果与分析通过对以上实验内容的实际操作和分析,得出以下结论:1. 在方程求解实验中,MATLAB给出的数值解与解析解基本吻合,验证了MATLAB在方程求解方面的高准确性和可靠性。

2. 在数值积分实验中,MATLAB给出的数值积分结果与解析积分结果基本吻合,验证了MATLAB在数值积分方面的高精度和稳定性。

3. 在常微分方程数值解实验中,MATLAB给出的数值解与解析解基本吻合,验证了MATLAB在常微分方程数值解方面的高准确性和可靠性。

五、结论与展望本实验通过MATLAB软件对数值分析方法进行了实验验证,得出了数值分析方法在不同数学问题上的高准确性和可靠性。

基于MATLAB数值分析实验报告

基于MATLAB数值分析实验报告

基于MATLAB数值分析实验报告班级:072115姓名:李凯学号:20111003943实验二:矩阵与向量运算实验目的:在MATLAB里,会对矩阵与向量进行加、减、数乘、求逆及矩阵特征值运算,以及矩阵的LU分解。

设A是一个n×n方阵,X是一个n维向量,乘积Y=AX可以看作是n维空间变换。

如果能够找到一个标量λ,使得存在一个非零向量X,满足:AX=λX (3.1)则可以认为线性变换T(X)=AX将X映射为λX,此时,称X 是对应于特征值λ的特征向量。

改写式(3.1)可以得到线性方程组的标准形式:(A-λI)X=0 (3.2)式(3.2)表示矩阵(A-λI)和非零向量X的乘积是零向量,式(3.2)有非零解的充分必要条件是矩阵(A-λI)是奇异的,即:det(A-λI)=0该行列式可以表示为如下形式:a11–λa12 (1)a21 a22 –λ…a2n =0 (3.3)…………A n1 a n2 …a nn将式(3.3)中的行列式展开后,可以得到一个n阶多项式,称为特征多项式:f(λ)=det(A-λI)=(-1)n(λn+c1λn-1+c2λn-2+…+c n-1λ+c n) (3.4) n阶多项式一共有n个根(可以有重根),将每个根λ带入式(3.2),可以得到一个非零解向量。

习题:求下列矩阵的特征多项式的系数和特征值λj:3 -1 0A= -1 2 -10-1 3解:在MATLAB中输入命令:A=【3 -1 0;-1 2 -1;0 -1 3】;c=poly(A)roots(c)得到:实验四:Lagrange插值多项式实验目的:理解Lagrange插值多项式的基本概念,熟悉Lagrange插值多项式的公式源代码,并能根据所给条件求出Lagrange插值多项式,理解龙格现象。

%功能:对一组数据做Lagrange插值%调用格式:yi=Lagran_(x,y,xi)%x,y:数组形式的数据表%xi:待计算y值的横坐标数组%yi:用Lagrange还擦之算出y值数组function fi=Lagran_(x,f,xi)fi=zeros(size(xi));np1=length(f);for i=1:np1z=ones(size(xi));for j=i:np1if i~=j,z=z.*(xi-x(j))/(x(i)-x(j));endendfi=fi+z*f(i);endreturn习题:已知4对数据(1.6,3.3),(2.7,1.22),(3.9,5.61),(5.6,2.94)。

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

曲线拟合的最小二乘法
1、
%采用二次多项式拟合
%a的输出为多项式各项系数
%b为拟合曲线各点函数值
%phi为输出的曲线拟合函数
x=0:5:55;
y=[0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.02 4.64]; a=polyfit(x,y,2)
b=polyval(a,x)
syms t
phi=a(1)*t^2+a(2)*t+a(3)
运行结果:
>> leastway
a =
-0.0024 0.2037 0.2305
b =
Columns 1 through 5
0.2305 1.1894 2.0293 2.7502 3.3521
Columns 6 through 10
3.8349
4.1987 4.4435 4.5693 4.5760
Columns 11 through 12
4.4637 4.2324
phi =
2
-0.00238051948051948162 t + 0.203690809190809258 t +
0.230467032967031749
2、
%采用三次多项式拟合
x=0:5:55;
y=[0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.02 4.64]; a=polyfit(x,y,3)
b=polyval(a,x)
syms t
phi=a(1)*t^3+a(2)*t^2+a(3)*t+a(4)
运行结果:
>> leastway
a =
0.0000 -0.0052 0.2634 0.0178
b =
Columns 1 through 5
0.0178 1.2087 2.1646 2.9113 3.4745
Columns 6 through 10
3.8800
4.1536 4.3211 4.4082 4.4407
Columns 11 through 12
4.4444 4.4450
phi =
3 2
0.0000343641543641541613 t - 0.00521556221556219567 t
+ 0.263398527398526872 t + 0.0178388278388323038
3、
%delta为拟合函数值与原函数值的误差
clc
clear
x=0:5:55;
y=[0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.02 4.64]; a=polyfit(x,y,3);
b=polyval(a,x);
for j=1:12
delta=b(j)-y(j)
end
syms t
phi=a(1)*t^3+a(2)*t^2+a(3)*t+a(4)
运行结果:
delta =
0.0178
delta =
-0.0613
delta =
0.0046
delta =
0.0513
delta =
0.0345
delta =
0.0100
delta =
0.0036
delta =
-0.0489
delta =
-0.1018
delta =
-0.1393
delta =
0.4244
delta =
-0.1950
phi =
3 2
0.0000343641543641541613 t - 0.00521556221556219567 t
+ 0.263398527398526872 t + 0.0178388278388323038
4、
%采用四次多项式拟合与三次多项式拟合进行比较
%输出其与原函数值的误差再与之前输出的三次与原函误差进行比较
%delta4为四次与原函的误差
clc
clear
x=0:5:55;
y=[0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.02 4.64]; a=polyfit(x,y,4);
b=polyval(a,x)
for j=1:12
delta4=b(j)-y(j)
end
phi=a(1)*t^4+a(2)*t^3+a(3)*t^2+a(4)*t+a(5)
运行结果:
b =
Columns 1 through 5
0.0604 1.1739 2.1220 2.8945 3.4900 Columns 6 through 10
3.9162
4.1898 4.3366 4.3914 4.3981 Columns 11 through 12
4.4095 4.4876
delta4 =
0.0604
delta4 =
-0.0961
delta4 =
-0.0380
delta4 =
0.0345
delta4 =
0.0500
0.0462
delta4 =
0.0398
delta4 =
-0.0334
delta4 =
-0.1186
delta4 =
-0.1819
delta4 =
0.3895
delta4 =
-0.1524
phi =
-6 4 3 0.602564102564107288 10 t - 0.0000319178969178975504t
2
- 0.00293227466977464190 t +0.238069314944314520t
+ 0.0604487179487194770
5、
%分别采用二次、四次、六次、九次多项式拟合并画出其曲线进行比较
clc
clear
x=0:5:55;
y=[0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.02 4.64]; a2=polyfit(x,y,2);
b2=polyval(a2,x);
a4=polyfit(x,y,4);
b4=polyval(a4,x);
a6=polyfit(x,y,6);
b6=polyval(a6,x);
a9=polyfit(x,y,9);
b9=polyval(a9,x);
plot(x,b2,'r-^');
hold on
plot(x,b4,'g->');
hold on
plot(x,b6,'b-*');
hold on
plot(x,b9,'y-+');
hold on
plot(x,y,'ko');
title('含碳量-时间曲线');
xlabel('时间t(分)');
ylabel('含碳量y(10^-^4)');
legend('二次拟合曲线','四次拟合曲线','六次拟合曲线','九次拟合曲线','原函数点',);
运行结果:
警告: 多项式未正确设置条件。

请添加具有不同X 值的点,减少多项式的阶数,或者尝试按照HELP
POLYFIT 中所述进行居中和缩放。

> In polyfit at 75
In leastway at 9
警告: 多项式未正确设置条件。

请添加具有不同X 值的点,减少多项式的阶数,或者尝试按照HELP
POLYFIT 中所述进行居中和缩放。

> In polyfit at 75
In leastway at 11
局部放大图。

相关文档
最新文档