完整word版matlab数值分析例题.docx
《完整word版matlab数值分析例题.docx》由会员分享，可在线阅读，更多相关《完整word版matlab数值分析例题.docx（12页珍藏版）》请在冰豆网上搜索。
完整word版matlab数值分析例题
1、在MATLAB中用Jacobi迭代法讨论线性方程组,
（1）给出Jacobi迭代法的迭代方程，并判定Jacobi迭代法求解此方程组是否收敛。
（2）若收敛，编程求解该线性方程组.
解
（1）：
A=[4-11;4—81;-215］%线性方程组系数矩阵
A=
4-11
4-81
—215
>>D=diag（diag（A））
D=
400
0—80
005
〉〉L=—tril（A,-1）％A的下三角矩阵
L=
000
—400
2—10
〉〉U=-triu（A，1）%A的上三角矩阵
U=
01—1
00—1
000
B=inv（D）*（L+U）%B为雅可比迭代矩阵
B=
00.2500—0。
2500
0.500000.1250
0。
4000—0.20000
〉〉r=eigs（B,1）%B的谱半径
r=
0。
3347〈1
Jacobi迭代法收敛。
（2）在matlab上编写程序如下:
A=[4—11;4-81;—215]；
〉〉b=［7—2115］'；
>〉x0=［000]’;
〉〉［x,k］=jacobi（A,b，x0，1e—7）
x=
2。
0000
4.0000
3。
0000
k=
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；%迭代次数
whilenorm（x-x0）>=eps
x0=x;
x=B＊x0+f；
k=k+1;
if（k〉=max1）
disp（’迭代超过300次，方程组可能不收敛’）；
return；
end
end
2、设有某实验数据如下：
序号
x
y
序号
x
y
1
—3
-3。
99
8
0.5
1。
3776
2
-2。
5
-3。
3011
9
1
1.5403
3
-2
-2。
4161
10
1。
5
1。
5707
4
—1。
5
—1。
4293
11
2
1.5839
5
—1
—0.4597
12
2。
5
1。
6989
6
—0.5
0。
37758
13
3
2.01
7
0
1
（1）在MATLAB中作图观察离散点的结构,用多项式拟合的方法拟合一个合适的多项式函数；
（2）在MATLAB中作出离散点和拟合曲线图。
解
（1）：
首先观察离散点的结构，matlab中的程序如下，
x=[—3-2.5-2-1。
5-1—0。
500.511.522。
53]；
>〉y=［—3.99-3。
3011-2。
4161-1.4293-0。
45970。
3775811.37761。
54031。
57071。
58391.69892。
01］；
>〉plot（x,y,'r＊’）
图形如下：
离散点近似如抛物线，所以用二次多项式拟合，所以matlab程序如下:
x=［-3-2.5—2-1.5-1—0。
500.511.522.53]；
〉〉y=[—3.99-3。
3011—2.4161—1。
4293—0。
45970。
3775811。
37761.54031.57071。
58391。
69892。
01]；
>〉s=polyfit（x,y，2）；
〉〉p=poly2str（s，'t’）
p=
-0。
22214t^2+1t+0.74384
（2）做出离散点与拟合曲线的程序如下：
x=[—3—2。
5—2-1。
5—1—0.500.511.522。
53];
〉>f=[—3.99—3。
3011—2。
4161—1.4293-0。
45970.3775811.37761。
54031.57071.58391。
69892。
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公式编程计算
（要求积分精度为
）
解：
matlab程序如下：
[I，step]=jfSimpson（'-x—x^2'，—1，2,2,1.0e—5）
I=
-4。
5000
step=
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
end
I=0；
switchtype
case1，
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;
case2,
n=2；
h=（b—a）/2；
I1=0;
I2=（subs（sym（f）,findsym（sym（f））,a）+subs（sym（f）,findsym（sym（f））,b））/h;
whileabs（I2—I1）>eps
n=n+1；
h=（b—a）/n;
I1=I2；
I2=0;
fori=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;
end
4、在MATLAB中用四阶Runge-Kutta法编程求解微分方程初值问题
，
并在MATLAB中画图比较方程的解析解与R—K解的结果.
解：
第1步：
首先先把经典RK4子程序编出来如下，用RK4。
m保存。
RK4。
mmatlab程序如下：
function［t,y］=RK4（func，t0，tt,y0，N,varagin）
%Rk方法计算一阶级微分方程组，
%由微分方程知识可以知道，对于高阶微分方程可以化为一阶微分方程组。
％初始时刻为t0，结束时为tt，初始时刻函数值为y0
%N为步数
ifnargin<4
N=100；
end
％如果没有输入N的话,那么默认计算步长为（tt-t0）/100
y（1，:
）=y0（：
）';
h=（tt—t0）/N；％步长
t=t0+[0:
N］’＊h;
fork=1:
N
f1=h＊feval（func，t（k），y（k，：
））；
f1=f1（：
）';
％RK中的k1
f2=h＊feval（func，t（k）+h/2,y（k，：
）+f1/2）；
f2=f2（:
）';
%RK中的k2
f3=h*feval（func,t（k）+h/2，y（k,:
）+f2/2）；
f3=f3（：
）’;
％RK中的k3
f4=h＊feval（func，t（k）+h,y（k,:
）+f3）;
f4=f4（：
）';
%RK中的k4
y（k+1,:
）=y（k，:
）+（f1+2*（f2+f3）+f4）/6；
％所求结果
End
第2步:
然后把微分方程的表达式用一个调试程序RK_fun表示，用RK_fun.m保存。
RK_fun.mmatlab程序如下：
％RK4方法的测试函数
functionf=RK_fun（t，y）
f=—y+t*t+3；
第3步：
把用RK方法计算的主函数写出，用RK_main.m保存。
RK_main.mmatlab程序如下：
%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
〉>yrk
yrk=
1.0000
1.5267
1.9647
2.3837
2.8352
3。
3575
3。
9789
4.7203
5。
5972
6.6213
7.8010
最后调用主函数得:
〉〉RK_main
得到分析结果如下图：
5、在MATLAB中利用牛顿切线迭代法计算非线性方程
在区间[1，2]上的一个根。
解，利用matlab写如下程序得:
root=newtonqxf（'x*x*x+4*x＊x-10'，1,2,1.0e—6）
root=
1.3652
附牛顿切线法程序如下：
functionroot=newtonqxf（f,a,b，eps）
%牛顿法求函数f在区间［a,b］上的一个零点
%f为函数名
％a为区间左端点
%b为区间右端点
％eps为根的精度
％root为求出的函数零点
if（nargin==3）
eps=1。
0e-6;
end
f1=subs（sym（f）,findsym（sym（f））,a）;
f2=subs（sym（f）,findsym（sym（f）），b）；
if（f1==0）
root=a;
end
if（f2==0）
root=b；
end
if（f1＊f2〉0）
disp（’两端点函数值乘积大于0！
’）；
return；
else
tol=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；
else
root=b-fb/dfb;
end
while（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
升级会员 