复制成功
  • 图案背景
  • 纯色背景

笔记

  • 2019-11-16
    为大人带来形象的羊生肖故事来历 为孩子带去快乐的生肖图画故事阅读
    谈谈怎样学好数学_苏步青-中学生文库

最新上海电力学院数值计算方法上机实习题

下载积分:2500

内容提示: 2017 数值计算方法 上机实习报告 学院: 专业: 班级: 姓名: 学号: 数值计算方法上机实习题 1. 设+=10 5dxxxInn, (1) 由递推公式nI In n151 +− =−,从0 =0.1823I , 1824 . 00= I 出发,计算20I ; (2) 20 =0I ,20 =10000I , 用nI In n51511 1+ − =− −,计算0I ; (3) 分析结果的可靠性及产生此现象的原因(重点分析原因)。 解:(1)分别令 I 0 的近似值为 0.1823、0.1824,MATLAB 程序如下: I=0.182...

文档格式:DOC| 浏览次数:36| 上传日期:2017-05-28 14:35:46| 文档星级:
2017 数值计算方法 上机实习报告 学院: 专业: 班级: 姓名: 学号: 数值计算方法上机实习题 1. 设+=10 5dxxxInn, (1) 由递推公式nI In n151 +− =−,从0 =0.1823I , 1824 . 00= I 出发,计算20I ; (2) 20 =0I ,20 =10000I , 用nI In n51511 1+ − =− −,计算0I ; (3) 分析结果的可靠性及产生此现象的原因(重点分析原因)。 解:(1)分别令 I 0 的近似值为 0.1823、0.1824,MATLAB 程序如下: I=0.1823; %题中的已知数据 for n=1:20; I=(-5)*I+1/n; %由递推公式所得 end fprintf('I20=%f\n',I) M=0.1824; %与I的计算结果形成对比 for i=1:20; M=(-5)*M+1/i; %由递推公式所得 end fprintf('M20=%f\n',M) %% 输出结果 I20=-2055816073.851284 M20=7480927090.212283 (2)分别令 I 20 的近似值为 0、10000,MATLAB 程序如下: I=0; %赋予I20的初始值 for n=0:19; I=(-1/5)*I+1/(5*(20-n)); %由递推公式所得 end fprintf('I0=%f\n',I) M=10000; for i=0:19; M=(-1/5)*M+1/(5*(20-i));%由递推公式所得 end fprintf('M0=%f\n',M) %% 输出结果 I0=0.182322 M0=0.182322 (3)分析: 由输出结果可看出第一种算法为不稳定算法,第二种算法为稳定算法。 由于误差 *0 0 0I I e − = 0 221 1*1*1 1*5 5 5 ) ( 5 )15 (15 e e e I InInI I I enn n n n n n n n n= = = − = + − − + − = − =− − − − − − 第一种算法为正向迭代算法,每计算一步误差增长 5 倍,虽然所给的初始值很接近,随着 n 的增大,误差也越来越大。 *n n nI I e − = n n n n n n n n ne I InInI e I I e51) (51)5151(5151* * *1 1 1= − = + − − + − = = − =− − − nne e )51(0= 第二种算法为倒向迭代算法,每计算一步误差缩小 5 倍,虽然所给的初始值之间差很多,随着 n 的增大,误差也越来越小。算法趋近稳定,收敛,可以选择这种算法。 2. 求方程 0 2 10 = − + x ex的近似根,要求4110 5−+× < −k kx x ,并比较计算量。 (1) 在[0,1]上用二分法; (2) 取初值 00= x ,并用迭代1021xkex−=+; (3) 加速迭代的结果; (4) 取初值 00= x ,并用牛顿迭代法; (5) 分析绝对误差。 解: (1)利用二分法,MATLAB 程序如下: %% 二分法程序 clear all clc a=0;b=1; f=@(x)(exp(x)+10*x-2); %@是定义函数句柄的运算符 c=(a+b)/2;%取区间中点 i=0;%分割次数 while abs(f(c))>5*10^(-4) %判断f(x)的精度是否满足要求 if f(a)*f(c)<0 b=c;c=(a+b)/2; elseif f(b)*f(c)<0 a=c;c=(b+a)/2; end i=i+1; end fprintf('二分法运算次数为%i\n',i) fprintf('二分法计算结果为%f\n',c) %% 输出结果 二分法运算次数为13 二分法计算结果为0.090515 (2)用题目中给出的迭代法,取初始值 x(1)=0,并用迭代 x(i)=(2-exp(x(i-1)))/10,MATLAB程序如下: %% 不动点迭代 clear all clc x0=0; x=x0; for k=1:10000 %规定迭代次数上限 y=(2-exp(x))/10; %迭代结果存到y中 if abs(x-y)<5*10^(-4) fprintf('初始值为x0%i\n迭代次数为%i\n',x0,k); break end x=y; if k==10000; fprintf('迭代次数超出上限%i\n',k); end end fprintf('迭代法计算结果为%f\n',y); %% 输出结果 初始值 x0 为 0 迭代次数为 4 迭代法计算结果为 0.090513 (3)利用加速迭代法,MATLAB 程序如下: %% 加速迭代法 x=0; f=@(x)(( 2-exp(x))./10); y=f(x); y0=f(y); y1=x-((y-x)^2)/(y0-2*y+x); i=1; while abs(y1-x)> 5*10^(-4) x=y1; y=f(x); y0=f(y); y1=x-((y-x)^2)/(y0-2*y+x); i=i+1; end %% 输出结果 y1= 0.090525 i=2 (4)牛顿迭代法: 取初始值 x=0,MATLAB 程序如下: %% 清空环境变量 clear all clc %% 牛顿迭代法 x=0; y=x-(exp(x)+10*x-2)./(exp(x)+10); i=1; while abs(x-y)>0.0005 x=y; y=x-(exp(x)+10*x-2)./(exp(x)+10); i=i+1; end %% 输出结果 y= 0.090525 i=2 (5)分析绝对误差 根据方程解求得对应的绝对误差如下表所示: 迭代方式 迭代次数 迭代结果 绝对误差 二分法迭代 11 0.0905762 5.1099e-05 不动点迭代 4 0.0905126 -1.2501e-05 加速迭代 2 0.0905251 -1.3073e-09 牛顿迭代 2 0.0905251 -1.3073e-09 通过上述表格可以看出,二分法运算了11次,不动点迭代方法运算了4次, Atiken加速迭代法运算了2次,牛顿迭代法运算2次.比较绝对误差可以发现 Atiken加速迭代和牛顿迭代的计算结果的绝对误差较小。下面就其原因进行分析: 我们知道1 2 1 2( ) ( ) g x g x L x x − ≤ − ,其中L为压缩常数,并且 0 1 L ≤ < 。 进行误差估计:1 0*1kkLx x x xL− ≤ −−,当L较小时,收敛较快,反之,当L很靠近1时,收敛很慢。若 1 L ≥ 时,则迭代不收敛。 由 ( ) 1 g x L ′ ≤ < ,( )1( )lim = ( 0, 1)!pkkke g xC C pe p∗+→∞= ≠ ≥ 常数 ,收敛阶为p。 二分法和不动点迭代为线性收敛;Atiken 迭代和牛顿迭代是平方收敛,是超线性收敛的。其中,二分法 L=0.5;用迭代1021xkex−=+,2( ) = = 0.2710 10 10x xe e eg x′  −′ ≤ ≈  ,比 0.5 小,因此收敛比二分法快。不动点迭代为线性收敛,Atiken 迭代速度与牛顿迭代速度都是超线性收敛的,因而收敛速度较快。相比之下,只有二分法的收敛速度较慢,绝对误差最大。 3.钢水包使用次数多以后,钢包的容积增大,数据如下: x 2 3 4 5 6 7 8 9 y 6.42 8.2 9.58 9.5 9.7 10 9.93 9.99 10 11 12 13 14 15 16 10.49 10.59 10.60 10.8 10.6 10.9 10.76试从中找出使用次数和容积之间的关系,计算均方差。(用ax byc x+=+拟合) 解:拟合曲线的模型是ax byc x+=+,将原模型变为 ) ( ) ( b ax x c y + = + ,采用非线性最小二程 法。按照最小二乘原理,应选取参数a,b,c使得表达式=+ − + =151 k2)) ( ) ( ( Z b ax x c yk k k达到极小值。具体的方法就是对Z关于a、b、c求偏导数,并置这些偏导数等于0,相对应的方程组如下所示: = + − + − =∂∂= + − + − =∂∂= + − + − =∂∂−−−0 )) ( ) (( 20 )) ( ) (( 20 )) ( ) (( 2151151151kk k k kkk k kkk k k kb ax y x c ycZb ax y x cbZb ax y x c xaZ 通过对上述的偏导数方程组进行整理,可以将之写成的 b AX = 的形式,利用 MATLA 充当计算器的作用,可以进行 X 的求解,对应的就是 a、b、c 三个参数的值。下面的程序就是进行求解方程组、画拟合曲线以及求取均方差的程序,具体如下: %求拟合方程,画拟合曲线 x=[2 3 4 5 6 7 8 9 10 11 12 13 14 15 16]; y=[6.42 8.2 9.58 9.5 9.7 10 9.93 9.99 10.49 10.59 10.60 10.8 10.6 10.9 10.76]; A=[sum(x.^2),sum(x),-sum(x.*y); sum(x),15,-sum(y); sum(x.*y),sum(y),-sum(y.^2)]; B=[sum(x.^2.*y);sum(x.*y);sum(y.^2.*x)]; f=A\B; a=f(1);b=f(2);c=f(3); z=linspace(2,50,48); Y=(f(1)*z+f(2))./(f(3)+z); plot(x,y,'r*',z,Y,'b-'); fprintf('a=%f\nb=%f\nc=%f\n',a,b,c); fprintf('拟合出来的方程式为:\n(%7.4f+x)y=%7.4f+%7.4fx\n',f(3),f(2),f(1)); %求均方差 for i=1:15 y1(i)=(f(2)+f(1)*x(i))/(f(3)+x(i)); end c=0; for i=1:15 c=c+(y(i)-y1(i))^2; end jfc=sqrt(c/15); fprintf('均方差为%f\n',jfc) %%运行结果: a=11.340048 b=-12.495325 c=-0.340291 拟合出来的方程式为: (-0.3403+x)y=-12.4953+11.3400x 均方差为 0.220812 对应的拟合曲线如下图: 分析总结:此方法利用线性方法求解非线性问题,避免了非线性误差较大的问题。 4.设− −− − −− − −− − −− − −− −=4 1 0 1 0 01 4 1 0 1 00 1 4 1 0 11 0 1 4 1 00 1 0 1 4 10 0 1 0 1 4A ,−−=625250b , b x = A 分析下列迭代法的收敛性,并求42110 −+≤ −k kx x 的近似解及相应的迭代次数。 (1) JACOBI 迭代; (2) GAUSS-SEIDEL 迭代; (3) SOR 迭代(取 0.1:0.1:1.9 ω = ,找到迭代步数最少的 * ω )。 解: (1)JACOBI 迭代 %% JACOBI函数 function tx=jacobi(A,b,imax,x0,tol) %初始值x0,次数imax,精度tol del=10^(-10); tx=[x0]; n=length(x0); for i=1:n dg=A(i,i); if abs(dg)<del disp('对角元太小');%防止出现溢出现象 return end end for k=1:imax %jacobi循环 for i=1:n sm=b(i); for j=1:n if j~=i sm=sm-A(i,j)*x0(j); end end x(i)=sm/A(i,i); %x(1)到x(n)即完成一次迭代 end tx=[tx;x];%矩阵中又加一行 if norm(x-x0)<tol return else x0=x; end end %% 主函数程序 clear all clc A=[4 -1 0 -1 0 0;-1 4 -1 0 -1 0;0 -1 4 -1 0 -1;-1 0 -1 4 -1 0;0 -1 0 -1 4 -1;0 0 -1 0 -1 4]; b=[0 5 -2 5 -2 6]; x0=[0 0 0 0 0 0]; imax=100;%迭代次数 tol=10^(-4);%精度 tx=jacobi(A,b,imax,x0,tol); for j=1:length(tx) fprintf('%d%f%f%f%f%f%f\n',j,tx(j,1),tx(j,2),tx(j,3),tx(j,4),tx(j,5),tx(j,6)); end %% 运行结果 10.0000000.0000000.0000000.0000000.0000000.000000 20.0000001.250000-0.5000001.250000-0.5000001.500000 30.6250001.0000000.5000001.0000000.5000001.250000 40.5000001.6562500.3125001.6562500.3125001.750000 50.8281251.5312500.7656251.5312500.7656251.656250 60.7656251.8398440.6796881.8398440.6796881.882813 70.9199221.7812500.8906251.7812500.8906251.839844 80.8906251.9252930.8505861.9252930.8505861.945313 90.9626461.8979490.9489751.8979490.9489751.925293 100.9489751.9651490.9302981.9651490.9302981.974487 110.9825741.9523930.9761961.9523930.9761961.965149 120.9761961.9837420.9674841.9837420.9674841.988098 130.9918711.9777910.9888951.9777910.9888951.983742 140.9888951.9924150.9848311.9924150.9848311.994448 150.9962081.9896390.9948201.9896390.9948201.992415 160.9948201.9964620.9929231.9964620.9929231.997410 170.9982311.9951670.9975831.9951670.9975831.996462 180.9975831.9983490.9966991.9983490.9966991.998792 190.9991751.9977450.9988731.9977450.9988731.998349 200.9988731.9992300.9984601.9992300.9984601.999436 210.9996151.9989480.9994741.9989480.9994741.999230 220.9994741.9996410.9992821.9996410.9992821.999737 230.9998201.9995090.9997551.9995090.9997551.999641 240.9997551.9998320.9996651.9998320.9996651.999877 250.9999161.9997710.9998861.9997710.9998861.999832 260.9998861.9999220.9998441.9999220.9998441.999943 270.9999611.9998930.9999471.9998930.9999471.999922 280.9999471.9999640.9999271.9999640.9999271.999973 290.9999821.9999500.9999751.9999500.9999751.999964 (2)GAUSS-SEIDEL 迭代 %% GAUSS-SEIDEL函数 function tx=gseidel(A,b,imax,x0,tol) %初始值x0,次数imax,精度tol del=10^-10; tx=[x0];n=length(x0); %tx是个二维矩阵,存储着每一步迭代的结果 for i=1:n dg=A(i,i); if abs(dg)<del disp('对角元太小'); return end end for k=1:imax %G-S循环 x=x0; for i=1:n sm=b(i); for j=1:n if j~=i sm=sm-A(i,j)*x(j); end end x(i)=sm/A(i,i); end tx=[tx;x]; if norm(x-x0)<tol return else x0=x; end end %% 主函数程序 clear all clc A=[4 -1 0 -1 0 0;-1 4 -1 0 -1 0;0 -1 4 -1 0 -1;-1 0 -1 4 -1 0;0 -1 0 -1 4 -1;0 0 -1 0 -1 4]; b=[0 5 -2 5 -2 6]; x0=[0 0 0 0 0 0]; imax=100;%迭代次数 tol=10^(-4);%精度 tx=gseidel(A,b,imax,x0,tol); for j=1:length(tx) fprintf('%d %f %f %f %f %f %f\n',j,tx(j,1),tx(j,2),tx(j,3),tx(j,4),tx(j,5),tx(j,6)) end %% 运行结果 1 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 2 0.000000 1.250000 -0.187500 1.203125 0.113281 1.481445 3 0.613281 1.384766 0.517334 1.560974 0.606796 1.781033 4 0.736435 1.715141 0.764287 1.776880 0.818263 1.895638 5 0.873005 1.863889 0.884102 1.893843 0.913342 1.949361 6 0.939433 1.934219 0.944356 1.949283 0.958216 1.975643 7 0.970875 1.968362 0.973322 1.975603 0.979902 1.988306 8 0.985991 1.984804 0.987178 1.988268 0.990344 1.994381 9 0.993268 1.992698 0.993837 1.994362 0.995360 1.997299 10 0.996765 1.996490 0.997038 1.997291 0.997770 1.998702 11 0.998445 1.998313 0.998577 1.998698 0.998928 1.999376 12 0.999253 1.999189 0.999316 1.999374 0.999485 1.999700 13 0.999641 1.999610 0.999671 1.999699 0.999752 1.999856 14 0.999827 1.999813 0.999842 1.999855 0.999881 1.999931 15 0.999917 1.999910 0.999924 1.999931 0.999943 1.999967 16 0.999960 1.999957 0.999964 1.999967 0.999973 1.999984 (3)SOR 迭代 %% SOR函数 function tx=sor(A,b,imax,x0,tol,w) %初始值x0,次数imax,精度tol del=10^-10; tx=[x0]; n=length(x0); for i=1:n dg=A(i,i); if abs(dg)<del disp('对角元太小'); return end end for k=1:imax %SOR循环 x=x0; for i=1:n sm=b(i); for j=1:n if j~=i sm=sm-A(i,j)*x(j); %先高斯迭代 end end x(i)=sm/A(i,i); x(i)=w*x(i)+(1-w)*x0(i); %比较高斯与超松弛迭代公式,补上省缺的项 end tx=[tx;x]; if norm(x-x0)<tol return else x0=x; end end %% 主函数程序 clear all clc A=[4 -1 0 -1 0 0;-1 4 -1 0 -1 0;0 -1 4 -1 0 -1;-1 0 -1 4 -1 0;0 -1 0 -1 4 -1;0 0 -1 0 -1 4]; b=[0 5 -2 5 -2 6]; x0=[0 0 0 0 0 0]; imax=100;迭代次数 tol=10^(-4);%精度 w=1.2;%给定不同的w,w=0.1:0.1:1.9,找出使SOR迭代法收敛速度最快 tx=sor(A,b,imax,x0,tol,w); for j=1:size(tx,1) fprintf('%d %f%f%f%f%f%f\n',j,tx(j,1),tx(j,2),tx(j,3),tx(j,4),tx(j,5),tx(j,6)) end %% 运行结果 1 0.0000000.0000000.0000000.0000000.0000000.000000 2 0.0000001.500000-0.1500001.4550000.2865001.840950 3 0.8865001.5069000.8708551.8221560.8937021.961177 4 0.8214171.9744120.9531531.9360500.9827511.988536 5 1.0088551.9885450.9833092.0052650.9981531.996732 6 0.9963721.9956411.0026291.9980940.9975092.000695 7 0.9988462.0005670.9992811.9990721.0005991.999825 8 1.0001231.9998870.9997792.0003360.9998951.999937 9 1.0000421.9999371.0001071.9999460.9999672.000035 10 0.9999572.0000220.9999791.9999821.0000181.999992 11 1.0000101.9999980.9999962.0000110.9999971.999999 W 值的范围在 0~2 之间,SOR 迭代才会收敛,令 w=0.1:0.1:1.9,找出 w*,使得 SOR迭代的收敛速度最快。W=0.1 时需要迭代 101 次;w=1.9 时需要迭代 101 次;w=1 时需要迭代 16 次;w=0.9 时需要迭代 20 次;w=1.1 时需要迭代 13 次;w=1.3 时需要迭代 13 次;w=1.2时需要迭代 11 次。由这几次试代可得出 w=1.2 时迭代次数最少,SOR 迭代法收敛速度最快。 总结: 由于 A 为严格对角占有矩阵,根据定理可知雅克比迭代和 GS 迭代都收敛,SOR 迭代收敛的必要条件是 0<w<2,即 SOR 迭代也收敛。三种迭代算法的结果分析如下: (1)JACOBI 迭代是根据 A 分解成上三角、下三角和对角阵,将线性方程组的求解归结为对角方程组求解过程的重复,思路简单易于编程。迭代次数比较多,收敛速度慢; (2)GAUSS-SEIDEL 迭代,是在 JACOBI 收敛的基础上将计算出来的未知量马上投入下一个迭代方程中使用,使得迭代的速度加快,效果更好。但 GS 迭代与 JACOBI 迭代的收敛域并不能相互包容,所以不能相互代替。但如果两者都收敛时,一般说 GS 迭代比 JACOBI迭代的收敛速度快。 (3)SOR 迭代为了加快收敛速度,在 GS 的基础上加入一修正参数即松弛因子 w,使得收敛速度更快。但选着不同的 w,收敛速度不一样,为了达到最好的效果,要选着最合适的松弛因子。 5.用逆幂迭代法求=1 1 11 2 31 3 6A 最接近于 11 的特征值和特征向量,准确到310 − 。 解: %% 逆幂迭代法 A=[6 3 1;3 2 1;1 1 1]; v(:,1)=[1;0;1]; tol=1e-3; p=11; L(1)=max(v(:,1)); u(:,1)=v(:,1)/L(1); k=1; B=(A-p*eye(3)); while k>0 v(:,k+1)=B\u(:,k); L(k+1)=max(v(:,k+1)); u(:,k+1)=v(:,k+1)/L(k+1); err(k+1)=abs(1/L(k+1)-1/L(k)); if err(k+1)<tol break else k=k+1; end end k u(:,k) 1/L(k)+p %% 运行结果 迭代次数 k =10 特征向量 = 4.4361 2.4363 1.0000 特征值 =7.8724 通过运行的结果可以看出,迭代运算了 10 次,最终求出的接近 11 的近似特征值为7.8724,对应的近似特征向量为[4.4361 2.4363 1.0000]’。 6.用经典 R-K 方法求解初值问题 (1)− + − = ′+ + − = ′x x y y yx y y ysin 2 cos 2 2sin 2 22 1 22 1 1, ] 10 , 0 [ ∈ x , ==3 ) 0 (2 ) 0 (21yy; (2)− + − = ′+ + − = ′x x y y yx y y ysin 999 cos 999 999 998sin 2 22 1 22 1 1, ] 10 , 0 [ ∈ x , ==3 ) 0 (2 ) 0 (21yy。 和精确解+ =+ =−−x e x yx e x yxxcos 2 ) (sin 2 ) (21比较,进行误差分析得到结论,图形显示精确解和数值解。 解: (1)MATLAB 程序如下: %% 用经典R-K方法求解初值问题 clc;clear; f=@(x,y1,y2)(-2*y1+y2+2*sin(x)) g=@(x,y1,y2) (y1-2*y2+2*cos(x)-2*sin(x)); h=0.1; y1(1)=2;y2(1)=3;x(1)=0; for i=1:100; K1=f(x(i),y1(i),y2(i)); L1=g(x(i),y1(i),y2(i)); K2=f(x(i)+0.5*h,y1(i)+0.5*h*K1,y2(i)+0.5*h*L1); L2=g(x(i)+0.5*h,y1(i)+0.5*h*K1,y2(i)+0.5*h*L1); K3=f(x(i)+0.5*h,y1(i)+0.5*h*K2,y2(i)+0.5*h*L2); L3=g(x(i)+0.5*h,y1(i)+0.5*h*K2,y2(i)+0.5*h*L2); K4=f(x(i)+h,y1(i)+h*K3,y2(i)+h*L3); L4=g(x(i)+h,y1(i)+h*K3,y2(i)+h*L3); x(i+1)=x(i)+h; y1(i+1)=y1(i)+h*(1/6)*(K1+2*K2+2*K3+K4); y2(i+1)=y2(i)+h*(1/6)*(L1+2*L2+2*L3+L4); end x=0:0.1:10; for i=1:101; Y1(i)=2*exp(-x(i))+sin(x(i)); Y2(i)=2*exp(-x(i))+cos(x(i)); end c=y1-Y1; d=y2-Y2; subplot(2,1,1) plot(x,y1,'r*',x,Y1,'b-','LineWidth',2) C=max(abs(y1-Y1)) C1=max(abs(y2-Y2)) legend('y1龙格库塔近似解','Y1精确解'); title('第一问曲线1精确解和近似解的对比'); ylabel('y1曲线 Y1曲线') subplot(2,1,2) plot(x,y2,'r*',x,Y2,'b-','LineWidth',2) legend('y2龙格库塔近似解','Y2精确解'); title('第一问曲线2精确解和近似解的对比'); ylabel('y2曲线 Y2曲线') %% 图形结果 误差结果: C = 7.3449e-06 C1 = 7.4367e-06 (2)MATLAB 程序如下: %% 用经典R-K方法求解初值问题 clc;clear; f=@(x,y1,y2)(-2*y1+y2+2*sin(x)); g=@(x,y1,y2)(998*y1-999*y2+999*cos(x)-999*sin(x)); h=0.001; y1(1)=2;y2(1)=3;x(1)=0; for i=1:10000; K1=f(x(i),y1(i),y2(i)); L1=g(x(i),y1(i),y2(i)); K2=f(x(i)+0.5*h,y1(i)+0.5*h*K1,y2(i)+0.5*h*L1); L2=g(x(i)+0.5*h,y1(i)+0.5*h*K1,y2(i)+0.5*h*L1); K3=f(x(i)+0.5*h,y1(i)+0.5*h*K2,y2(i)+0.5*h*L2); L3=g(x(i)+0.5*h,y1(i)+0.5*h*K2,y2(i)+0.5*h*L2); K4=f(x(i)+h,y1(i)+h*K3,y2(i)+h*L3); L4=g(x(i)+h,y1(i)+h*K3,y2(i)+h*L3); x(i+1)=x(i)+h; y1(i+1)=y1(i)+h*(1/6)*(K1+2*K2+2*K3+K4); y2(i+1)=y2(i)+h*(1/6)*(L1+2*L2+2*L3+L4); end x=0:0.001:10; Y1=2*exp(-x)+sin(x); Y2=2*exp(-x)+cos(x);%¾«È·½â subplot(211) plot(x,y1,'r*',x,Y1,'b-','LineWidth',2) legend('y1龙格库塔近似解','Y1精确解'); title('第二问曲线1精确解和近似解的对比'); ylabel('y1曲线 Y1曲线') subplot(2,1,2) plot(x,y2,'r*',x,Y2,'b-','LineWidth',2) legend('y2龙格库塔近似解','Y2精确解'); title('第二问曲线2精确解和近似解的对比'); ylabel('y2曲线 Y2曲线') C=max(abs(y1-Y1)) C1=max(abs(y2-Y2)) %% 图形结果 误差结果: C = 2.4263e-11 C1 = 2.3547e-08 (3)结论 四阶龙格库塔法所求出的解的精确度很高,很接近精确解。证明该方法是一种有效的算法。考虑到方程组的刚性问题,为了保证算法的稳定性,必须将步长限制在较小的范围内。第一组数据的刚性比为 3 大于 1,为病态方程组,要求其步长限制在较小的范围内。第二组数据的刚性比为 1000 也大于 1,为病态方程组,故步长也要限制在较小范围内。 7.用有限差分法求解边值问题(h=0.1),并图形显示。 = = −= + − ′ ′1 ) 1 ( ) 1 (0 ) 1 (2y yy x y. 解: %% 有限差分法求解边值问题 %% 追赶算法 function x=chase(a,b,c,f) %求解线性方程组Ax=f,其中A使三对角阵 %a是矩阵A的小对角线元素a(1)=0 %b是矩阵A的对角线元素 %c是矩阵A的上对角线元素 %f是方程组的右端向量 N=length(f); x=zeros(1,N);y=zeros(1,N); d=zeros(1,N);u= zeros(1,N); %预处理 d(1)=b(1); for i=1:N-1 u(i)=c(i)/d(i); d(i+1)=b(i+1)-a(i+1)*u(i); end %追的过程 y(1)=f(1)/d(1); for i=2:N y(i)=(f(i)-a(i)*y(i-1))/d(i); end %赶的过程 x(N)=y(N); for i=N-1:-1:1 x(i)=y(i)-u(i)*x(i+1); end %% 主程序 h=0.1; N=(1-(-1))/h; x(1)=-1; x(N-1)=1; y(1)=1; y(N-1)=1; for i=1:N-1 X(i)=x(1)+i*h; q(i)=1+(X(i)^2); B(i)=2/(h^2)+q(i); end for i=1:N-2 C(i)=-1/(h^2); end g(1)=0+1/(h^2); g(N-1)=0+1/(h^2); for i=2:N-2 g(i)=0; end y=chase([0,C],B,[C,0],g) y1=[1,y,1]; X=-1:h:1; plot(X,y1); title('有限差分法'); xlabel('输入值x'); ylabel('y曲线'); %% 输出结果及图形 y = Columns 1 through 10 0.9096 0.8357 0.7754 0.7267 0.6879 0.6577 0.6352 0.6195 0.6103 0.6073 Columns 11 through 19 0.6103 0.6195 0.6352 0.6577 0.6879 0.7267 0.7754 0.8357 0.9096 总结:步长分别取万分之一和十分之一时对比结果如上图,采用万分之一步长时可以近似为精确解,采用十分之一步长时可以认为是近似解,步长越小,运算越多。

关注我们

关注微信公众号

您选择了以下内容