数值计算方法源程序.docx

1、数值计算方法源程序计算物理实验目 录1.迭代法求解方程 12.主序列消元法求解方程组 23.追赶法求解方程组 34.牛顿插值 45.复合梯形公式求积分 56.复合Simpson公式求解积分 67. 欧拉法求初值问题 68. 预估校正法求初值问题 89. 四阶龙格库塔法求初值问题 910.备注： 101.迭代法求解方程（1）M文件编辑为：function diedfun(x0,wc) %定义m文件函数 x0是初始值，wc是误差k=0;x=x0+2*wc;w=1;while wwc&k x0=2; %输入初始值 wc=0.0001; %输入误差 diedfun(x0,wc) %调用函数计算结果是:

2、x(0)=2.828498 x(1)=2.942277 x(2)=2.980697 x(3)=2.993559 x(4)=2.997852 x(5)=2.999284 x(6)=2.999761 x(7)=2.999920 x(8)=2.9999731即，迭代9次就达到我们所要的精度。所得的解是约等于3.0（2）M文件编辑：function diedfun(x0,wc) %定义m文件函数k=0;x=x0+2*wc;w=1;while wwc&k x0=2; %输入初始值 wc=0.0001; %输入允许误差 diedfun(x0,wc) %调用函数运行结果如下：x(49996)=-0.9936

3、66x(49997)=-1.006314x(49998)=-0.993667x(49999)=-1.006313迭代50000次后仍不能达到我们所预设的精度,收敛性很差。 2.主序列消元法求解方程组M文件编辑：gxy.m function gxy(a,b) %定义主序列消元函数m=length(b);x=zeros(1,m);for k=1:m %计算每列的最大值，并进行换行运算，保证在下三角形中对角占优 for g=1:m c=a(g,g); for l=g+1:m if c a=1,1,1;12,-3,3;-18,3,-1; %定义系数，用矩阵表示 b=6,15,-15; %定义右侧值，用

4、一维矩阵表示 gxy(a,b) %求解方程组计算结果是：方程组的计算结果是：x(1)=1.000000 x(2)=2.000000 x(3)=3.0000003.追赶法求解方程组M文件编辑：zhuigan.mfunction zhuigan(a,b,c,f) %定义m文件函数n=length(f);for i=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;endx(n)=f(n)/b(n); %赶的过程（也就是回代的过程）for i=n-1:-1:1 x(i)=(f(i)-c(i)*x(i+1)/b

5、(i);endfor i=1:n fprintf(x(%d)=%.6fn,i,x(i);endcommand windows 中输入： 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.000000x(2)=1.000000x(3)=1.000000x(4)=1.0000004.牛顿插值文件编辑：newtonin.mfunction newtonin(x,y,xx) %定义m文件函数m=lengt

6、h(x);n=length(y);r=length(xx);p=zeros(1,m);if m=n error(cuowu!);endfor i=2:m %计算各插值系数 k=n; while k=i y(k)=(y(k-1)-y(k)/(x(k-i+1)-x(k); k=k-1; end endfor j=1:rp(j)=y(n) ;i=n-1;while i=1 %计算插值结果 p(j)=p(j)*(xx(j)-x(i)+y(i) ; i=i-1;endfprintf(n%f的结果是%f,xx(j),p(j);end在command windows中输入： x=0,0.3,0.6,0.9,

7、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.9665880.500000的结果是0.9588490.750000的结果是0.9088541.000000的结果是0.8414675.复合梯形公式求积分M文件编辑：（1）.tixing.mfunction tixing(a,b,wc) %定义m文件程序h=b-a;l=1;s1=(tifun(a)+ti

8、fun(b)*h;w=1;while wwc %利用精度控制积分的精度 s=0; for i=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(积分的结果是：%fn,s);（2）tifun.mfunction t=tifun(x) %m文件函数定义的被积函数 t=1/(1+x*x);command windows中输入： a=0; %积分下限 b=1; %积分上限 wc=0.00001; %精度 tixing(a,b,wc); %调用积分函数计算

9、结果是：积分的结果是：0.785389分析：要达到一定精度需要运行的时间比较长。6.复合Simpson公式求解积分M文件编辑：（1）.simpson.mfunction simpson(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;while wwc %利用精度控制积分过程 f3=0; for i=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=

10、s; n=n*2; h=h/2; f2=f2+f3;endfprintf(积分计算的结果是：%f,s);（2）simfun.mfunction t=simfun(x) %m文件定义的被积函数 t=1/(1+x*x); %被积函数的返回值command windows 中输入： a=0; %积分下限 b=1; %积分上限 wc=0.00001; %精度 simpson(a,b,wc); %调用积分函数运行结果是：积分计算的结果是：0.785398计算的时间比梯形的要短的多。7. 欧拉法求初值问题M文件编辑：（1）.ola.mfunction ola(a,b,n,y0) %定义m文件的积分函数h=

11、(b-a)/n;n=n+1;x=zeros(1,n);y=x;y(1)=y0; %给y赋初值for i=2:n %此循环求各x点的值和对应的y的值 x(i)=x(i-1)+h; y(i)=y(i-1)+h*olafun(x(i-1),y(i-1); endplot(x,y,-b*);hold on; %画模拟的曲线q=exp(x);plot(x,q);hold on; %画指数函数曲线（2）olafun.mfunction d=olafun(x,y); %m文件定义的f(x)d=x*y;command windows 中输入： a=0; %（a，b）设定被求的区间 b=1; n=10; %设置

12、分割次数 y0=1; %设置初值 ola(a,b,n,y0) %调用函数运行结果如图所示：8. 预估校正法求初值问题M文件编辑：（1）.yuce.mfunction yuce(a,b,m,y0) %m文件定义函数h=(b-a)/mm=m+1;t=zeros(1,m);t(1)=a;for i=2:m t(i)=t(i-1)+h;endy=zeros(1,m);y(1)=y0;for i=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)

13、+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;endfor i=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;endplot(t,y,-b

14、*); hold on; %画模拟的曲线q=exp(t);plot(t,q);hold on; %画指数函数的曲线（2）yf.mfunction t=yf(x,y) %m文件定义函数t=x*y;command windows 中输入： a=0; %（a,b）定义计算的区间 b=1; n=10; %分割次数 y0=1; %定义初值 yuce(a,b,n,y0) %调用函数运行结果如图所示：9. 四阶龙格库塔法求初值问题M文件编辑：（1）.rk.mfunction rk(a,b,n,y0) %m文件定义函数h=(b-a)/n;n=n+1;t=zeros(1,n);y=t;x(1)=a;for i=

15、2:n t(i)=t(i-1)+h;endy(1)=y0;for i=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;endplot(t,y,-b*);hold on; %绘制模拟曲线q=exp(t);plot(t,q);hold on; %绘制指数函数曲线（2）rkf.mfunction z=rkf(x

16、,y) %m文件定义函数z=x*y;command windows 中输入： a=0; %（a,b）定义的计算区间 b=1; n=10; %分割次数 y0=1; %定义初值 rk(a,b,n,y0) %调用函数运行结果如图所示：10.备注：（1）所有的程序均是用MATLAB编写，并已运行通过。（2）简要使用说明：打开MATLAB，出现命令窗口（command windows）在命令窗口中输入要编辑的m文件，如第九题，则输入：edit rk,然后回车，出现m文件编辑器，在其中输入编写的函数程序，其中放入函数名称一定要和文件名称一致，并且可以带参数，如第九题：函数名为rk(a,b,n,y0)。然后点击保存，切记：一定要保存，才会有用，不然修改的结果将得不到体现各循环，条件，判断，均以end作为结尾标记。在命令窗口中输入各初始值，然后回车，就会得到我们所要的结果。%后面的是注释部分，不会运行的。