forq=1:
m
d=a（l,q）;
a（l,q）=a（g,q）;
a（g,q）=d;
end
d=b（l）;
b（l）=b（g）;
b（g）=d;
end
end
end
fori=k:
m
c=a（i,k）;
forj=k:
m
a（i,j）=a（i,j）/c;%将每行的对角元素化为1
end
b（i）=b（i）/c;
end
fori=k:
m
forj=k:
m-1
a（j+1,i）=a（j+1,i）-a（k,i）;%将下三角化为0
end
end
fori=k:
m-1
b（i+1）=b（i+1）-b（k）;
end
end
x（m）=b（m）;
fori=m-1:
-1:
1%回代求解过程
x（i）=b（i）;
forj=i+1:
m
x（i）=x（i）-a（i,j）*x（j）;
end
x（i）=x（i）/a（i,i）;
end
fprintf（'方程组的计算结果是：
\n'）
fori=1:
m
fprintf（'x（%d）=%.6f\t',i,x（i））;
end
commandwindows中输入：
>>a=[1,1,1;12,-3,3;-18,3,-1];%定义系数，用矩阵表示
>>b=[6,15,-15];%定义右侧值，用一维矩阵表示
>>gxy（a,b）%求解方程组
计算结果是：
方程组的计算结果是：
x
（1）=1.000000x
（2）=2.000000x（3）=3.000000
3.追赶法求解方程组
M文件编辑：
zhuigan.m
functionzhuigan（a,b,c,f）%定义m文件函数
n=length（f）;
fori=2:
n%追的过程
t=a（i）/b（i-1）;
a（i）=0;
b（i）=b（i）-c（i-1）*t
f（i）=f（i）-f（i-1）*t;
end
x（n）=f（n）/b（n）;%赶的过程（也就是回代的过程）
fori=n-1:
-1:
1
x（i）=（f（i）-c（i）*x（i+1））/b（i）;
end
fori=1:
n
fprintf（'x（%d）=%.6f\n',i,x（i））;
end
commandwindows中输入：
>>a=[0,-1,-1,-1];%定义第一列函数
>>b=[2,2,2,2];%定义第二列函数
>>c=[-1,-1,-1,0];%定义第三列函数
>>f=[1,0,0,1];%定义右侧系数
>>zhuigan（a,b,c,f）%调用函数解方程
出现运行结果为：
x
（1）=1.000000
x
（2）=1.000000
x（3）=1.000000
x（4）=1.000000
4.牛顿插值
Ｍ文件编辑：
newtonin.m
functionnewtonin（x,y,xx）%定义m文件函数
m=length（x）;
n=length（y）;
r=length（xx）;
p=zeros（1,m）;
ifm~=n
error（'cuowu!
'）;
end
fori=2:
m%计算各插值系数
k=n;
whilek>=i
y（k）=（y（k-1）-y（k））/（x（k-i+1）-x（k））;
k=k-1;
end
end
forj=1:
r
p（j）=y（n）;
i=n-1;
whilei>=1%计算插值结果
p（j）=p（j）*（xx（j）-x（i））+y（i）;
i=i-1;
end
fprintf（'\n%f的结果是%f',xx（j）,p（j））;
end
在commandwindows中输入：
>>x=[0,0.3,0.6,0.9,1.2];%x数组
>>y=[1.000006,0.9850674,0.9410708,0.8703632,0.7766992];%对应于x的y值
>>xx=[0.45,0.5,0.75,1];%要求的一系列的值
>>newtonin（x,y,xx）%调用牛顿插值公式
运行结果如下：
0.450000的结果是0.966588
0.500000的结果是0.958849
0.750000的结果是0.908854
1.000000的结果是0.841467
5.复合梯形公式求积分
M文件编辑：
（1）.tixing.m
functiontixing（a,b,wc）%定义m文件程序
h=b-a;
l=1;s1=（tifun（a）+tifun（b））*h;
w=1;
whilew>wc%利用精度控制积分的精度
s=0;
fori=1:
l
x=a+（i-1/2）*h;
s=s+tifun（（x））;
end
s=h*s/2;
s=s+s1/2;
h=h/2;l=l*2;
w=abs（s-s1）;
s1=s;
end
fprintf（'积分的结果是：
%f\n',s）;
（2）tifun.m
functiont=tifun（x）%m文件函数定义的被积函数
t=1/（1+x*x）;
commandwindows中输入：
>>a=0;%积分下限
>>b=1;%积分上限
>>wc=0.00001;%精度
>>tixing（a,b,wc）;%调用积分函数
计算结果是：
积分的结果是：
0.785389
分析：
要达到一定精度需要运行的时间比较长。
6.复合Simpson公式求解积分
M文件编辑：
（1）.simpson.m
functionsimpson（a,b,wc）%定义的m文件积分函数
h=b-a;
f1=simfun（a）+simfun（b）;
f2=simfun（（b+a）/2）;
s0=（f1+4*f2）*h/6;
h=h/4;
w=1;
n=2;
whilew>wc%利用精度控制积分过程
f3=0;
fori=1:
n
f3=f3+simfun（a+（2*i-1）*h）;
end
s=（f1+2*f2+4*f3）*h/3;
w=abs（s0-s）*15;
s0=s;n=n*2;h=h/2;f2=f2+f3;
end
fprintf（'积分计算的结果是：
%f',s）;
（2）simfun.m
functiont=simfun（x）%m文件定义的被积函数
t=1/（1+x*x）;%被积函数的返回值
commandwindows中输入：
>>a=0;%积分下限
>>b=1;%积分上限
>>wc=0.00001;%精度
>>simpson（a,b,wc）;%调用积分函数
运行结果是：
积分计算的结果是：
0.785398
计算的时间比梯形的要短的多。
7.欧拉法求初值问题
M文件编辑：
（1）.ola.m
functionola（a,b,n,y0）%定义m文件的积分函数
h=（b-a）/n;
n=n+1;
x=zeros（1,n）;
y=x;
y
（1）=y0;%给y赋初值
fori=2:
n%此循环求各x点的值和对应的y的值
x（i）=x（i-1）+h;
y（i）=y（i-1）+h*olafun（x（i-1）,y（i-1））;
end
plot（x,y,'--b*'）;holdon;%画模拟的曲线
q=exp（x）;
plot（x,q）;holdon;%画指数函数曲线
（2）olafun.m
functiond=olafun（x,y）;%m文件定义的f（x）
d=x*y;
commandwindows中输入：
>>a=0;%（a，b）设定被求的区间
>>b=1;
>>n=10;%设置分割次数
>>y0=1;%设置初值
>>ola（a,b,n,y0）%调用函数
运行结果如图所示：
8.预估—校正法求初值问题
M文件编辑：
（1）.yuce.m
functionyuce（a,b,m,y0）%m文件定义函数
h=（b-a）/m
m=m+1;
t=zeros（1,m）;
t
（1）=a;
fori=2:
m
t（i）=t（i-1）+h;
end
y=zeros（1,m）;
y
（1）=y0;
fori=2:
4%此循环求解从2到4各x值所对应的y值
k1=h*yf（t（i-1）,y（i-1））;
k2=h*yf（t（i-1）+h/2,y（i-1）+k1/2）;
k3=h*yf（t（i-1）+h/2,y（i-1）+k2/2）;
k4=h*yf（t（i-1）+h,y（i-1）+k3）;
y（i）=y（i-1）+（k1+2*k2+2*k3+k4）/6;
end
fori=5:
m%此循环求解剩下的x值所对应的y值
y（i）=y（i-1）+h*（55*yf（t（i-1）,y（i-1））-59*yf（t（i-2）,y（i-2））+37*yf（t（i-3）,y（i-3））-9*yf（t（i-4）,y（i-4）））/24;
y（i）=y（i-1）+h*（9*yf（t（i）,y（i））+19*yf（t（i-1）,y（i-1））-5*yf（t（i-2）,y（i-2））+yf（t（i-3）,y（i-3）））/24;
end
plot（t,y,'--b*'）;holdon;%画模拟的曲线
q=exp（t）;
plot（t,q）;holdon;%画指数函数的曲线
（2）yf.m
functiont=yf（x,y）%m文件定义函数
t=x*y;
commandwindows中输入：
>>a=0;%（a,b）定义计算的区间
>>b=1;
>>n=10;%分割次数
>>y0=1;%定义初值
>>yuce（a,b,n,y0）%调用函数
运行结果如图所示：
9.四阶龙格—库塔法求初值问题
M文件编辑：
（1）.rk.m
functionrk（a,b,n,y0）%m文件定义函数
h=（b-a）/n;
n=n+1;
t=zeros（1,n）;
y=t;
x
（1）=a;
fori=2:
n
t（i）=t（i-1）+h;
end
y
（1）=y0;
fori=2:
n%利用初值，计算各点的值
k1=rkf（t（i-1）,y（i-1））;
k2=rkf（t（i-1）+h/2,y（i-1）+h*k1/2）;
k3=rkf（t（i-1）+h/2,y（i-1）+h*k2/2）;
k4=rkf（t（i-1）+h,y（i-1）+h*k3）;
y（i）=y（i-1）+（k1+2*k2+2*k3+k4）*h/6;
end
plot（t,y,'--b*'）;holdon;%绘制模拟曲线
q=exp（t）;
plot（t,q）;holdon;%绘制指数函数曲线
（2）rkf.m
functionz=rkf（x,y）%m文件定义函数
z=x*y;
commandwindows中输入：
>>a=0;%（a,b）定义的计算区间
>>b=1;
>>n=10;%分割次数
>>y0=1;%定义初值
>>rk（a,b,n,y0）%调用函数
运行结果如图所示：
10.备注：
（1）所有的程序均是用MATLAB编写，并已运行通过。
（2）简要使用说明：
①打开MATLAB，出现命令窗口（commandwindows）
②在命令窗口中输入要编辑的m文件，如第九题，则输入：
editrk,然后回车，出现m文件编辑器，在其中输入编写的函数程序，其中放入函数名称一定要和文件名称一致，并且可以带参数，如第九题：
函数名为rk（a,b,n,y0）。
然后点击保存，切记：
一定要保存，才会有用，不然修改的结果将得不到体现
③各循环，条件，判断，均以end作为结尾标记。
④在命令窗口中输入各初始值，然后回车，就会得到我们所要的结果。
⑤%后面的是注释部分，不会运行的。
升级会员 