计算方法实验报告Word文件下载.docx
- 文档编号:20056376
- 上传时间:2023-01-16
- 格式:DOCX
- 页数:18
- 大小:167.63KB
计算方法实验报告Word文件下载.docx
《计算方法实验报告Word文件下载.docx》由会员分享,可在线阅读,更多相关《计算方法实验报告Word文件下载.docx(18页珍藏版)》请在冰豆网上搜索。
%计算A的维数为m
x(m,2)=0;
%初始化一个m行2列的矩阵,记录相邻两次迭代结果
x(:
2)=x0;
fori=[1:
1:
N],%迭代最高次数为N
ifmax(abs(x(:
1)-x(:
2)))<
E,%判断是否满足要求的精度
break;
else
x(:
1)=x(:
2);
forj=[1:
n],
ifj==1,
s1=0;
fore1=[2:
m],
s1=s1+A(1,e1)*x(e1,1);
%求和
end
x(j,2)=(b(1,1)-s1)/A(1,1);
elseifj==m,
s2=0;
fore2=[1:
m-1],
s2=s2+A(m,e2)*x(e2,2);
x(j,2)=(b(m,1)-s2)/A(m,m);
s3=0;
fore3=[1:
j-1],
s3=s3+A(j,e3)*x(e3,2);
fore4=[j+1:
s3=s3+A(j,e4)*x(e4,1);
x(j,2)=(b(j,1)-s3)/A(j,j);
end
end
x1=x(:
【结果分析与讨论】
在matlab命令窗口中输入:
>
A=[5-1-1-1;
-110-1-1;
-1-15-1;
-1-1-110;
];
b=[-4;
12;
8;
34;
x0=[1;
1;
1];
GS(A,b,x0,5,0.000001)
ans=
0.996552514803221
1.998599188612620
2.998420412137192
3.999357211555304
评价:
相对于简单迭代法,赛德尔迭代法收敛速度很快,此程序可以对任意维的方程组求指定的精度的迭代结果,实用性比较强。
实验二最小二乘法拟合数据
熟悉最小二乘法拟合数据的方法
1.了解matlab语言
2.给出数据:
-1.00
0.50
-0.50
-0.25
0.00
10.2209
0.3295
0.8826
1.4392
2.0003
0.25
0.75
1.00
2.5645
3.1334
3.7061
4.2836
用一次、二次多项式利用最小二乘法拟合这些数据,试写出正规方程组,并求出最小平方逼近多项式。
最小二乘法拟合算法:
输入数据列向量
与
,拟合多项式的次数
计算
的列数为
(其中
表示矩阵
的第1行第1列,以下的
等类似);
利用
分解法求解
,得到
次最小平方逼近多项式的系数与多项式;
S6:
在同一坐标系中绘出原数据点与次
多项式逼近曲线。
实验程序:
functiony=fmintwo(a,m,x)%调用拟合多项式的函数
y=a(1,1);
fore1=[1:
y=y+(a(e1+1,1))*(x.^e1);
functiona=mintwo(x,y,m)%最小二乘法拟合程序
[n,p]=size(x);
a(m+1,1)=0;
A(m+1,m+1)=0;
b(m+1,1)=0;
c(2*m+1,1)=0;
c(1,1)=n;
fork=[2:
2*m+1],
s1=0;
fore1=[1:
s1=s1+(x(e1,1))^(k-1);
end
c(k,1)=s1;
m+1],
forj=[1:
A(i,j)=c(i+j-1,1);
%求正规方程组的系数矩阵
A
forh=[1:
s2=0;
fore2=[1:
s2=s2+y(e2,1)*((x(e2,1))^(h-1));
b(h,1)=s2;
%求正规方程组的
向量
b
[LDa]=ja2911(A,b);
%调用
分解法求解正规矩阵的程序
z=[min(x):
0.1:
max(x)];
plot(x,y,'
*'
z,fmintwo(a,m,z),'
r'
);
%绘出原数据表示的点和拟合曲线
附
分解法求系数矩阵为对称矩阵的方程组程序:
function[LDx]=ja2911(A,b)
[n,n]=size(A);
L=eye(size(A));
D=zeros(size(A));
fori=[2:
A(i,1)=A(i,1)/A(1,1);
forj=[2:
ifi>
j;
s1=s1+A(e1,e1)*A(i,e1)*A(j,e1);
A(i,j)=(A(i,j)-s1)/A(j,j);
elseifi==j,
i-1],
s2=s2+A(e2,e2)*A(i,e2)*A(i,e2);
A(i,j)=A(i,j)-s2;
D(h,h)=A(h,h);
fort=[1:
k-1],
L(k,t)=A(k,t);
y(n,1)=0;
x(n,1)=0;
y(1,1)=b(1,1)/D(1,1);
fora=[2:
s3=0;
fore3=[1:
a-1],
s3=s3+L(a,e3)*D(e3,e3)*y(e3,1);
y(a,1)=(b(a,1)-s3)/D(a,a);
x(n,1)=y(n,1);
forbb=[n-1:
-1:
1],
s4=0;
fore4=[bb+1:
s4=s4+L(e4,bb)*x(e4,1);
x(bb,1)=y(bb,1)-s4;
xx=[-1;
-0.75;
-0.5;
-0.25;
0;
0.25;
0.5;
0.75;
yy=[-0.2209;
0.3295;
0.8826;
1.4392;
2.0003;
2.5645;
3.1334;
…
3.7061;
4.2836;
mintwo(xx,yy,1)
A=
9.0000000000000000
03.750000000000000
b=
18.118299999999998
8.443674999999999
2.013144444444444
2.251646666666666
输出的图像:
mintwo(xx,yy,2)
9.00000000000000003.750000000000000
03.7500000000000000
3.75000000000000002.765625000000000
7.586956250000000
2.000100432900433
.0313********
一次多项式拟合的正规方程组为:
一次拟合多项式为:
二次多项式拟合的正规方程组为:
二次拟合多项式为:
此程序实现了对数据利用最小二乘法进行
次多项式逼近,并且可以输出原数据表示的点和拟合曲线在同一坐标系中的图像,可以很直观的进行观察拟合的程度。
实验三龙贝格求积分
熟悉龙贝格求积分的方法
2.用龙贝格方法计算积分:
要求误差不超过
。
龙贝格求积分算法:
S1:
输入区间左右端点
和精度
S7:
,则输出
,并且结束;
时,当
则转到s5;
则无输出,并且结束。
functiony=fromberg(x)%可被调用的积分函数
y=2*exp(-(x^2))/((pi)^(1/2));
functionI=romberg(a,b,E)%龙贝格求积分程序
T(50,1)=0;
C(50,1)=0;
R(50,1)=0;
h(50,1)=0;
h(1,1)=(b-1)/2;
T(1,1)=h(1,1)*(fromberg(a)+fromberg(b));
fork1=[2:
5],
h(k1,1)=(b-a)/(2^(k1-1));
g1=0;
2^(k1-2)],
g1=g1+fromberg(a+(2*e1-1)*h(k1,1));
T(k1,1)=(T(k1-1,1)/2)+h(k1,1)*g1;
%计算
fori1=[1:
4],
S(i1,1)=(4*T(i1+1,1)-T(i1,1))/3;
fori2=[1:
3],
C(i2,1)=(16*S(i2+1,1)-S(i2,1))/15;
fori3=[1:
2],
R(i3,1)=(64*C(i3+1,1)-C(i3,1))/63;
fork=[1:
40],
ifabs(R(k+1,1)-R(k,1))<
=E,%误差判断
I=R(k+1,1);
break;
else%如果不符合精度要求,则自动将区间分半继续计算
h(k+5,1)=(b-a)/(2^(k+4));
g2=0;
2^(k+3)],
g2=g2+fromberg(a+(2*e2-1)*h(k+5,1));
T(k+5,1)=(T(k+4,1)/2)+h(k+5,1)*g2;
S(k+4,1)=(4*T(k+5,1)-T(k+4,1))/3;
C(k+3,1)=(16*S(k+4,1)-S(k+3,1))/15;
R(k+2,1)=(64*C(k+3,1)-C(k+2,1))/63;
romberg(0,1,10^(-5))
0.842693582047394
龙贝格求积分方法精确度比较高,该程序很好的实现了龙贝格方法,但是程序的某些地方还可以简化,使程序精简。
实验四标准四阶R-K方法
熟悉标准四阶R-K的方法
2.用标准四阶R-K求解初值问题
的解在-0.9,-0.8,-0.7的近似值(按四位小数计算)。
、
的初值
、初始步长
、最大迭代次数
、允许误差
,则将
,输出
,将
,并转向s4;
,无输出,并且结束。
functionfk=fRK(x,y)%微分方程的
函数
fk=(x^2)-(y^2);
functionx=RK(a,b,ya,h)%用R-K方法求微分方程数值解的程序
[m,n]=size([a:
h:
b]);
x(n,2)=0;
1)=[a:
b]'
;
x(1,2)=ya;
k(4,1)=0;
k(1,1)=h*fRK(x(i-1,1),x(i-1,2));
k(2,1)=h*fRK(x(i-1,1)+h/2,x(i-1,2)+k(1,1)/2);
k(3,1)=h*fRK(x(i-1,1)+h/2,x(i-1,2)+k(2,1)/2);
k(4,1)=h*fRK(x(i-1,1)+h,x(i-1,2)+k(3,1));
x(i,2)=x(i-1,2)+(k(1,1)+2*k(2,1)+2*k(3,1)+k(4,1))/6;
function[xh1]=WRK(a,b,ya,h,N,E)%R-K的步长自动选择程序
m(N,1)=0;
x1(N,1)=0;
N],
[lm(i,1)]=size([a:
h/(2^(i-1)):
end;
y=RK(a,b,ya,h);
x1(1,1)=y(m(1,1),2);
y=RK(a,b,ya,h/2);
x1(2,1)=y(m(2,1),2);
forj=[3:
ifabs(x1(j-1,1)-x1(j-2,1))<
=E,
x=RK(a,b,ya,h/(2^(j-2)));
h1=h/(2^(j-2));
else
y=RK(a,b,ya,h/(2^(j-1)));
x(j,1)=y(m(3,1),2);
(该程序输出的是满足精度要求的区间分点的数值解和及相应步长)
[xh]=WRK(-1,-0.5,0,0.1,20,0.00001)
x=
-1.0000000000000000
-0.9500000000000000.047503044462813
-0.9000000000000000.090047822422416
-0.8500000000000000.127736329909335
-0.8000000000000000.160727765024206
-0.7500000000000000.189********6077
-0.7000000000000000.213483856644815
-0.6500000000000000.233766*********
-0.6000000000000000.250369794944545
-0.5500000000000000.263601751827377
-0.5000000000000000.273776792527207
h=
0.050000000000000
所以微分方程的解在-0.9,-0.8,-0.7的近似值分别为0.09004、0.16072、0.21348
标准四阶R-K方法的精确度比较高,而且算法容易实现。
次程序实现了标准四阶R-K算法,通过多次调用函数的形式使程序得到了精简。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 计算方法 实验 报告