修改稿 雅克比.docx
《修改稿 雅克比.docx》由会员分享，可在线阅读，更多相关《修改稿 雅克比.docx（13页珍藏版）》请在冰豆网上搜索。
修改稿雅克比
第一题牛顿插值法
function[A,C,L]=newtoncz（X,Y）
n=length（X）
A=zeros（n,n）
A（:
1）=Y';
forj=2:
n
fori=j:
n
A（i,j）=（A（i,j-1）-A（i-1,j-1））/（X（i）-X（i-j+1））;
end
end
C=A（n,n）;
fork=（n-1）:
-1:
1
C=conv（C,poly（X（k）））;
d=length（C）;
C（d）=C（d）+A（k,k）;
end
L=poly2sym（C）;
end
在matlab中输入
X=[0.20.40.60.81.0];
Y=[0.980.920.810.640.38];
[A,C,L]=newtoncz（X,Y）
>>M=polyval（C,0.75）
第二题最小二乘法
>>clear
>>x=[00.91.93.03.95.0];
>>y=[010305080110];
>>a=polyfit（x,y,2）
a=
2.248811.0814-0.5834
>>x1=[0:
0.5:
5.0];
>>y1=a（3）+a
（2）*x1+a
（1）*x1.^2;
>>plot（x,y,'-r'）;
>>holdon
>>plot（x1,y1,'*'）
>>xi=[00.91.93.03.95.0];
>>y=[1.752.453.814.807.008.60];
>>n=length（xi）;
>>x=[0:
0.5:
5.0];
f=-0.5834+11.0814*xi+2.2488*xi.^2;
>>fy=abs（f-y）
fy=
2.33348.761424.779448.100069.8383102.4436
>>fy2=fy.^2;
>>E3=sum（fy2）/n
E3=
3.0637e+03
第三题文件复合梯形
1fht
functionT_n=fht（a,b,n）
h=（b-a）/n;
fork=0:
n
x（k+1）=a+k*h;
ifx（k+1）==0
x（k+1）=10^（-10）;
end
end
T_1=h/2*（fx（x
（1））+fx（x（n+1）））;
fori=2:
n
F（i）=h*fx（x（i））;
end
T_2=sum（F）
T_n=T_1+T_2
文件2fx.m
functiony=fx（x）
y=sin（x）./x;
end
matlab中输入的
>>fht（0,1,8）
T_2=
0.8306
T_n=
0.9457
ans=
0.9457
>>formatlong
>>fht（0,1,8）
T_2=
0.830598927032208
T_n=
0.945690863582701
ans=
0.945690863582701
>>formatshort
>>fht（0,1,8）
复合S_P_S.m文件
functionS_n=S_P_S（a,b,n）
h=（b-a）/n;
fork=0:
n
x（k+1）=a+k*h;
x_k（k+1）=x（k+1）+1/2*h;
if（x（k+1）==0）|（x_k（k+1）==0）
x（k+1）=10^（-10）;
x_k（k+1）=10^（-10）;
end
end
S_1=（fx（x
（1））+fx（x（n+1）））*（h/6）;
fori=2:
n
F_1（i）=h/3*fx（x（i））;
end
forj=1:
n
F_2（j）=2*h/3*fx（x_k（j））;
end
S_2=sum（F_1）+sum（F_2）;
S_n=S_1+S_2;
文件名fx.x
functiony=fx（x）
y=sin（x）./x;
end
在matlab输入
>>S_1=s_p_s（0,1,8）
S_1=
0.9461
第四题用牛顿法解方程（223）：
用牛顿法求解方程
在1.5附近的根（要求误差
）。
functiony=newton（x0）;
x1=x0-fc（x0）/dfc（x0）;
n=1;
while（abs（x1-x0）>=1.0e-6）&（n<=100000000）
x0=x1;
x1=x0-fc（x0）/dfc（x0）;
n=n+1;
end
x1
n
end
第二个文件
functiony=fc（x）
y=x.^3-x.^2-1
functiony=dfc（x）
y=3*x.^2-2*x
在matlab中输入
>>newton（1.5）
y=
0.1250
y=
3.7500
y=
0.0039
y=
3.5200
y=
4.0700e-06
y=
3.5126
y=
4.5599e-12
y=
3.5126
x1=
1.4656
n=
4
5、雅克比迭代法与高斯—赛德尔迭代法（187）：
function[x,i,G]=kb（A,b,max,eps）
D=diag（diag（A））;
ID=inv（D）;
J=ID*（-A+D）;
f=ID*b;
x0=zeros（rank（A）,1）;
x=x0;
G（1,:
）=x0';
fori=1:
max
k=x;
x=J*x+f;
G（i+1,:
）=x';
n=norm（（x-k）,inf）;
ifn<=eps
break;
end
end
G
i
J
x
在matlab中输入
A=[1023;2101;3110];
b=[141120]';
max=100;
eps=1.0e-006;
[x,i,G]=kb（A,b,max,eps）
G=
000
1.4000000000000001.1000000000000002.000000000000000
0.5800000000000000.6200000000000001.470000000000000
0.8350000000000000.8370000000000001.764000000000000
0.7034000000000000.7566000000000001.665800000000000
0.7489400000000000.7927400000000001.713320000000000
0.7274560000000000.7788800000000001.696044000000000
0.7354108000000000.7849044000000001.703875200000000
0.7318565600000000.7825303200000001.700886320000000
0.7332280400000000.7835400560000001.702190000000000
0.7326349888000000.7831353920000001.701677582400000
0.7328696468800000.7833052440000001.701895964160000
0.7327701619520000.7832364742080001.701808581536000
0.7328101306976000.7832651094560001.701845303993600
0.7327933869107200.7832534434611201.701830449845120
0.7328001763542400.7832582776333441.701836639580672
0.7327973525991300.7832563007710851.701834119330393
0.7327985040466650.7832571175471351.701835164143153
0.7327980272476270.7832567827763521.701834737031287
i=
18
J=
0-0.200000000000000-0.300000000000000
-0.2000000000000000-0.100000000000000
-0.300000000000000-0.1000000000000000
x=
0.732798027247627
0.783256782776352
1.701834737031287
x=
0.732798027247627
0.783256782776352
1.701834737031287
i=
18
G=
000
1.4000000000000001.1000000000000002.000000000000000
0.5800000000000000.6200000000000001.470000000000000
0.8350000000000000.8370000000000001.764000000000000
0.7034000000000000.7566000000000001.665800000000000
0.7489400000000000.7927400000000001.713320000000000
0.7274560000000000.7788800000000001.696044000000000
0.7354108000000000.7849044000000001.703875200000000
0.7318565600000000.7825303200000001.700886320000000
0.7332280400000000.7835400560000001.702190000000000
0.7326349888000000.7831353920000001.701677582400000
0.7328696468800000.7833052440000001.701895964160000
0.7327701619520000.7832364742080001.701808581536000
0.7328101306976000.7832651094560001.701845303993600
0.7327933869107200.7832534434611201.701830449845120
0.7328001763542400.7832582776333441.701836639580672
0.7327973525991300.7832563007710851.701834119330393
0.7327985040466650.7832571175471351.701835164143153
0.7327980272476270.7832567827763521.701834737031287
六、常微分方程初值问题数值解法（318）：
用经典四阶K-R方法解初值问题
>>
步长分别取h=0.1,0.025,0.01，计算各步长下，数值解在x=0.4,0.6,0.8处的误差，指出误差值随步长的变化规律（已知真解
）。
function[]=LGKT1（h,x0,y0,X,Y,eps）
formatlong
h=input（'h='）;
x0=input（'x0='）;
y0=input（'y0='）;
disp（'输入的范围是：
'）;
X=input（'X='）;
Y=input（'Y='）;
disp（'输入的精度是：
'）;
eps=input（'eps='）;
i=1;
x1=0;
k1=0;
k2=0;
k3=0;
k4=0;
n=（Y-X）/h;
fori=1:
1:
n
x1=x0+h;
k1=-50*y0+50*x0.^2+2*x0;
k2=（-（x0+h/2）*（y0+h/2*k1）^2）;
k3=（-（x0+h/2）*（y0+h/2*k2）^2）;
k4=（-（x1）*（y0+h*k3）^2）;
y1=y0+h/6*（k1+2*k2+2*k3+k4）;
x0=x1;
y0=y1;
y=（1./3）*exp（-50*x0）+x0.^2;
fy=abs（y-y1）;
if（fy-eps）<=0
break
end
fprintf（'结果=%.3f,%.7f,%.7f,%.7f\n',x1,fy,y1,y）';
end
end
在MATLAB中输入
>>LGKT1
h=0.2
x0=0
y0=1/3
输入的范围是：
X=0
Y=1
输入的精度是：
eps=1.0e-8
结果=0.200,0.2754850,-0.2354698,0.0400151
结果=0.400,0.0467358,0.2067358,0.1600000
结果=0.600,0.2068270,0.1531730,0.3600000
结果=0.800,0.1825481,0.4574519,0.6400000
结果=1.000,0.3357547,0.6642453,1.0000000
升级会员 