微分方程数值解法课程设计许.doc

微分方程数值解法课程设计----A组
计算022班3号 许**
1. 用复化梯形计算积分
解：
算法思想：
根据复化梯形公式:取步长，
有，写出程序，如下所示：
N=100;
a=0;b=10;
h=(b-a)/N;
T=zeros(0,N);
for j=1:N
T(j)=0.5*h*[exp(-(a+h*(j-1)).^2)+exp(-(a+h*j).^2)]
t=sum(T)
end
运行结果：t = 0.8862
2. 用Euler、改进的Euler法、梯形法、R-K法解

解：
根据算法思想：
（1）Euler法

（2）改进的Euler法

（3）梯形法

（4）R-K法

我们可以写出程序如下所示：
format long;
N=100;
h=1/N;%步长
U=zeros(1,N+1);%euler
T=zeros(1,N+1);%真解
X=zeros(1,N+1);
R=zeros(4,N+1);%误差
G=zeros(1,N+1);%改进eule
O=zeros(1,N+1);%梯形法
K=zeros(1,N+1)%R-K法
c=zeros(1,4);
a=0;
for j=2:N+1;
X(j)=h*(j-1);
end
U(1)=1
T(1)=1;
G(1)=1;
X(1)=0;
O(1)=1;
K(1)=1;
for j=1:N
U(j+1)=U(j)+h*(U(j)+X(j)); %euler法公式
a=G(j)+h*(G(j)+X(j));
G(j+1)=G(j)+h/2*(G(j)+X(j)+a+X(j+1)); %改进的euler法公式
b=O(j)+h*(O(j)+X(j));
for i=1:10
O(j+1)=O(j)+h/2*(G(j)+X(j)+b+X(j+1)); %梯形法公式
b=O(j+1);
end
c(1)=K(j)+X(j);
c(2)=K(j)+h/2*c(1)+(X(j)+h/2);
c(3)=K(j)+h/2*c(2)+(X(j)+h/2);
c(4)=K(j)+h*c(3)+(X(j)+h);
K(j+1)=K(j)+h/6*(c(1)+2*c(2)+2*c(3)+c(4));%四阶R-K法公式
T(j+1)=-X(j+1)-1+2*exp(X(j+1)); %真实值
R(1,j+1)=U(j+1)-T(j+1); %EULER法误差
R(2,j+1)=G(j+1)-T(j+1); %改进的EULER法误差
R(3,j+1)=O(j+1)-T(j+1); %梯形法误差
R(4,j+1)=K(j+1)-T(j+1); %四阶R-K法误差
end
真实值T：
1 1.0101 1.0204 1.0309 1.0416 1.0525 1.0637 1.075 1.0866 1.0983 1.1103 1.1226 1.135 1.1477 1.1605 1.1737 1.187 1.2006 1.2144 1.2285 1.2428 1.2574 1.2722 1.2872 1.3025 1.3181 1.3339 1.3499 1.3663 1.3829 1.3997 1.4169 1.4343 1.4519 1.4699 1.4881 1.5067 1.5255 1.5446 1.564 1.5836 1.6036 1.6239 1.6445 1.6654 1.6866 1.7081 1.73 1.7521 1.7746 1.7974 1.8206 1.8441 1.8679 1.892 1.9165 1.9413 1.9665 1.9921 2.018 2.0442 2.0709 2.0979 2.1252 2.153 2.1811 2.2096 2.2385 2.2678 2.2974 2.3275 2.358 2.3889 2.4202 2.4519 2.484 2.5166 2.5495 2.5829 2.6168 2.6511 2.6858 2.721 2.7566 2.7927 2.8293 2.8663 2.9038 2.9418 2.9803 3.0192 3.0586 3.0986 3.139 3.18 3.2214 3.2634 3.3059 3.3489 3.3925 3.4366
欧拉法U值：
1 1.01 1.0202 1.0306 1.0412 1.052 1.063 1.0743 1.0857 1.0974 1.1092 1.1213 1.1337 1.1462 1.1589 1.1719 1.1852 1.1986 1.2123 1.2262 1.2404 1.2548 1.2694 1.2843 1.2995 1.3149 1.3305 1.3464 1.3626 1.379 1.3957 1.4127 1.4299 1.4474 1.4652 1.4832 1.5015 1.5202 1.5391 1.5582 1.5777 1.5975 1.6176 1.638 1.6586 1.6796 1.7009 1.7225 1.7445 1.7667 1.7893 1.8122 1.8354 1.8589 1.8828 1.907 1.9316 1.9565 1.9818 2.0074 2.0334 2.0597 2.0864 2.1135 2.1409 2.1687 2.1969 2.2255 2.2544 2.2838 2.3135 2.3437 2.3742 2.4051 2.4365 2.4683 2.5004 2.533 2.5661 2.5995 2.6334 2.6678 2.7025 2.7378 2.7734 2.8096 2.8462 2.8832 2.9208 2.9588 2.9973 3.0362 3.0757 3.1157 3.1561 3.1971 3.2385 3.2805 3.323 3.3661 3.4096
改进的欧拉法G值：
1 1.0101 1.0204 1.0309 1.0416 1.0525 1.0637 1.075 1.0866 1.0983 1.1103 1.1226 1.135 1.1477 1.1605 1.1737 1.187 1.2006 1.2144 1.2285 1.2428 1.2573 1.2721 1.2872 1.3025 1.318 1.3338 1.3499 1.3662 1.3828 1.3997 1.4168 1.4342 1.4519 1.4699 1.4881 1.5066 1.5255 1.5446 1.5639 1.5836 1.6036 1.6239 1.6445 1.6654 1.6866 1.7081 1.73 1.7521 1.7746 1.7974 1.8206 1.844 1.8678 1.892 1.9165 1.9413 1.9665 1.992 2.0179 2.0442 2.0708 2.0978 2.1252 2.1529 2.181 2.2095 2.2384 2.2677 2.2974 2.3275 2.3579 2.3888 2.4201 2.4518 2.4839 2.5165 2.5495 2.5829 2.6167 2.651 2.6858 2.7209 2.7566 2.7927 2.8292 2.8663 2.9038 2.9417 2.9802 3.0191 3.0586 3.0985 3.1389 3.1799 3.2213 3.2633 3.3058 3.3488 3.3924 3.4365
梯形法O值：
1 1.0101 1.0204 1.0309 1.0416 1.0525 1.0637 1.075 1.0866 1.0984 1.1103 1.1226 1.135 1.1477 1.1606 1.1737 1.187 1.2006 1.2144 1.2285 1.2428 1.2574 1.2722 1.2872 1.3025 1.3181 1.3339 1.3499 1.3663 1.3829 1.3997 1.4169 1.4343 1.4519 1.4699 1.4881 1.5067 1.5255 1.5446 1.564 1.5837 1.6036 1.6239 1.6445 1.6654 1.6866 1.7082 1.73 1.7522 1.7746 1.7975 1.8206 1.8441 1.8679 1.892 1.9165 1.9414 1.9665 1.9921 2.018 2.0442 2.0709 2.0979 2.1252 2.153 2.1811 2.2096 2.2385 2.2678 2.2974 2.3275 2.358 2.3889 2.4202 2.4519 2.484 2.5166 2.5495 2.583 2.6168 2.6511 2.6858 2.721 2.7567 2.7927 2.8293 2.8663 2.9038 2.9418 2.9803 3.0192 3.0587 3.0986 3.139 3.18 3.2214 3.2634 3.3059 3.3489 3.3925 3.4366
R-K法K值：
1 1.0101 1.0204 1.0309 1.0416 1.0525 1.0637 1.075 1.0866 1.0983 1.1103 1.1226 1.135 1.1477 1.1605 1.1737 1.187 1.2006 1.2144 1.2285 1.2428 1.2574 1.2722 1.2872 1.3025 1.3181 1.3339 1.3499 1.3663 1.3829 1.3997 1.4169 1.4343 1.4519 1.4699 1.4881 1.5067 1.5255 1.5446 1.564 1.5836 1.6036 1.6239 1.6445 1.6654 1.6866 1.7081 1.73 1.7521 1.7746 1.7974 1.8206 1.8441 1.8679 1.892 1.9165 1.9413 1.9665 1.9921 2.018 2.0442 2.0709 2.0979 2.1252 2.153 2.1811 2.2096 2.2385 2.2678 2.2974 2.3275 2.358 2.3889 2.4202 2.4519 2.484 2.5166 2.5495 2.5829 2.6168 2.6511 2.6858 2.721 2.7566 2.7927 2.8293 2.8663 2.9038 2.9418 2.9803 3.0192 3.0586 3.0986 3.139 3.18 3.2214 3.2634 3.3059 3.3489 3.3925 3.4366
误差R值：
第一组
0 -0.00010033 -0.00020268 -0.00030707 -0.00041353 -0.00052209 -0.00063279 -0.00074566 -0.00086072 -0.00097802 -0.0010976 -0.0012194 -0.0013436 -0.0014702 -0.0015992 -0.0017306 -0.0018645 -0.0020008 -0.0021398 -0.0022813 -0.0024254 -0.0025722 -0.0027217 -0.002874 -0.003029 -0.0031868 -0.0033475 -0.0035111 -0.0036777 -0.0038472 -0.0040198 -0.0041954 -0.0043742 -0.0045561 -0.0047412 -0.0049296 -0.0051213 -0.0053163 -0.0055147 -0.0057166 -0.0059219 -0.0061308 -0.0063433 -0.0065595 -0.0067793 -0.0070029 -0.0072303 -0.0074615 -0.0076966 -0.0079358 -0.0081789 -0.0084261 -0.0086775 -0.008933 -0.0091928 -0.0094569 -0.0097254 -0.0099983 -0.010276 -0.010558 -0.010844 -0.011135 -0.011431 -0.011732 -0.012038 -0.012349 -0.012664 -0.012985 -0.013311 -0.013642 -0.013979 -0.014321 -0.014668 -0.015021 -0.015379 -0.015743 -0.016113 -0.016489 -0.01687 -0.017258 -0.017651 -0.018051 -0.018457 -0.01887 -0.019288 -0.019714 -0.020146 -0.020584 -0.02103 -0.021482 -0.021941 -0.022407 -0.02288 -0.023361 -0.023849 -0.024344 -0.024847 -0.025358 -0.025876 -0.026402 -0.026936
第二组
0 -3.3417e-007 -6.7505e-007 -1.0228e-006 -1.3774e-006 -1.739e-006 -2.1078e-006 -2.4838e-006 -2.8672e-006 -3.258e-006 -3.6564e-006 -4.0624e-006 -4.4763e-006 -4.8981e-006 -5.3278e-006 -5.7658e-006 -6.212e-006 -6.6665e-006 -7.1296e-006 -7.6014e-006 -8.0818e-006 -8.5712e-006 -9.0696e-006 -9.5772e-006 -1.0094e-005 -1.062e-005 -1.1156e-005 -1.1702e-005 -1.2257e-005 -1.2822e-005 -1.3398e-005 -1.3983e-005 -1.458e-005 -1.5186e-005 -1.5804e-005 -1.6432e-005 -1.7071e-005 -1.7722e-005 -1.8384e-005 -1.9057e-005 -1.9742e-005 -2.0439e-005 -2.1148e-005 -2.1869e-005 -2.2603e-005 -2.3349e-005 -2.4108e-005 -2.4879e-005 -2.5664e-005 -2.6462e-005 -2.7273e-005 -2.8098e-005 -2.8937e-005 -2.979e-005 -3.0657e-005 -3.1539e-005 -3.2435e-005 -3.3346e-005 -3.4272e-005 -3.5213e-005 -3.617e-005 -3.7142e-005 -3.8131e-005 -3.9135e-005 -4.0156e-005 -4.1193e-005 -4.2247e-005 -4.3318e-005 -4.4407e-005 -4.5513e-005 -4.6636e-005 -4.7778e-005 -4.8938e-005 -5.0116e-005 -5.1313e-005 -5.2529e-005 -5.3765e-005 -5.502e-005 -5.6294e-005 -5.7589e-005 -5.8904e-005 -6.024e-005 -6.1596e-005 -6.2974e-005 -6.4373e-005 -6.5794e-005 -6.7237e-005 -6.8703e-005 -7.0191e-005 -7.1702e-005 -7.3236e-005 -7.4794e-005 -7.6376e-005 -7.7982e-005 -7.9613e-005 -8.1269e-005 -8.2949e-005 -8.4656e-005 -8.6388e-005 -8.8147e-005 -8.9932e-005
第三组
0 1.6834e-007 3.3755e-007 5.076e-007 6.7848e-007 8.5018e-007 1.0227e-006 1.196e-006 1.3701e-006 1.5449e-006 1.7205e-006 1.8968e-006 2.0739e-006 2.2516e-006 2.43e-006 2.6091e-006 2.7888e-006 2.9692e-006 3.1501e-006 3.3317e-006 3.5138e-006 3.6964e-006 3.8796e-006 4.0633e-006 4.2475e-006 4.4321e-006 4.6172e-006 4.8027e-006 4.9885e-006 5.1747e-006 5.3613e-006 5.5481e-006 5.7353e-006 5.9227e-006 6.1103e-006 6.2981e-006 6.486e-006 6.6741e-006 6.8623e-006 7.0506e-006 7.2389e-006 7.4272e-006 7.6155e-006 7.8037e-006 7.9918e-006 8.1798e-006 8.3676e-006 8.5552e-006 8.7425e-006 8.9295e-006 9.1162e-006 9.3025e-006 9.4884e-006 9.6738e-006 9.8587e-006 1.0043e-005 1.0227e-005 1.041e-005 1.0592e-005 1.0774e-005 1.0955e-005 1.1135e-005 1.1314e-005 1.1492e-005 1.1669e-005 1.1846e-005 1.2021e-005 1.2194e-005 1.2367e-005 1.2538e-005 1.2708e-005 1.2877e-005 1.3044e-005 1.3209e-005 1.3373e-005 1.3535e-005 1.3696e-005 1.3854e-005 1.4011e-005 1.4166e-005 1.4319e-005 1.4469e-005 1.4618e-005 1.4764e-005 1.4908e-005 1.5049e-005 1.5188e-005 1.5324e-005 1.5458e-005 1.5588e-005 1.5716e-005 1.5841e-005 1.5963e-005 1.6082e-005 1.6198e-005 1.631e-005 1.6419e-005 1.6524e-005 1.6626e-005 1.6724e-005 1.6818e-005
第四组
0 -1.6693e-012 -3.3724e-012 -5.1097e-012 -6.8812e-012 -8.6882e-012 -1.053e-011 -1.2409e-011 -1.4324e-011 -1.6277e-011 -1.8267e-011 -2.0295e-011 -2.2363e-011 -2.447e-011 -2.6617e-011 -2.8805e-011 -3.1034e-011 -3.3305e-011 -3.5619e-011 -3.7976e-011 -4.0376e-011 -4.2821e-011 -4.5311e-011 -4.7847e-011 -5.0429e-011 -5.3058e-011 -5.5735e-011 -5.846e-011 -6.1235e-011 -6.4059e-011 -6.6934e-011 -6.986e-011 -7.2838e-011 -7.5869e-011 -7.8954e-011 -8.2093e-011 -8.5287e-011 -8.8537e-011 -9.1844e-011 -9.5208e-011 -9.8631e-011 -1.0211e-010 -1.0565e-010 -1.0926e-010 -1.1292e-010 -1.1665e-010 -1.2044e-010 -1.2429e-010 -1.2821e-010 -1.322e-010 -1.3625e-010 -1.4038e-010 -1.4457e-010 -1.4883e-010 -1.5316e-010 -1.5756e-010 -1.6204e-010 -1.6659e-010 -1.7122e-010 -1.7592e-010 -1.807e-010 -1.8556e-010 -1.905e-010 -1.9551e-010 -2.0061e-010 -2.058e-010 -2.1106e-010 -2.1641e-010 -2.2185e-010 -2.2738e-010 -2.3299e-010 -2.3869e-010 -2.4449e-010 -2.5037e-010 -2.5635e-010 -2.6243e-010 -2.686e-010 -2.7487e-010 -2.8124e-010 -2.8771e-010 -2.9428e-010 -3.0095e-010 -3.0773e-010 -3.1461e-010 -3.216e-010 -3.287e-010 -3.3591e-010 -3.4323e-010 -3.5067e-010 -3.5822e-010 -3.6588e-010 -3.7367e-010 -3.8157e-010 -3.8959e-010 -3.9774e-010 -4.0601e-010 -4.1441e-010 -4.2293e-010 -4.3159e-010 -4.4037e-010 -4.4929e-010
3. 用差分格式计算两点边值问题

解：
算法思想：根据有

%-u''(t)+u(t)=[(1/20*x^4)-6]*x;
%u(0)=0，U(1)=21/20,U(x)=1/20*x^5+x^3
N=100;
X=zeros(1,N+1);
F=zeros(N-1,1);
h=1/N;
A=zeros(N-1,N-1);
U=zeros(N-1,1);
R=zeros(N-1,1);
X=0:h:1;
T=zeros(N-1,1);
A(1,1)=2+h^2;
A(1,2)=-1;
A(N-1,N-1)=2+h^2;
A(N-1,N-2)=-1;
for i=2:N-2
A(i,i-1)=-1;
A(i,i)=2+h^2;
A(i,i+1)=-1;
end
for j=1:N-1
F(j)=h^2*[(1/20*X(j+1)^4)-6]*X(j+1);
T(j)=1/20*X(j+1)^5+X(j+1)^3;
end
U=inv(A)*F;
R=U-T
运行结果：

U值：
-0.0089337 -0.017862 -0.026781 -0.035684 -0.044566 -0.053424 -0.06225 -0.071041 -0.079791 -0.088494 -0.097147 -0.10574 -0.11428 -0.12275 -0.13114 -0.13946 -0.1477 -0.15585 -0.16391 -0.17187 -0.17973 -0.18748 -0.19511 -0.20263 -0.21002 -0.21729 -0.22442 -0.23141 -0.23826 -0.24495 -0.2515 -0.25788 -0.26409 -0.27013 -0.276 -0.28168 -0.28717 -0.29248 -0.29758 -0.30248 -0.30717 -0.31164 -0.31589 -0.31992 -0.32372 -0.32727 -0.33059 -0.33365 -0.33646 -0.33902 -0.3413 -0.34331 -0.34505 -0.3465 -0.34766 -0.34853 -0.3491 -0.34936 -0.34931 -0.34894 -0.34824 -0.34721 -0.34585 -0.34415 -0.34209 -0.33968 -0.33691 -0.33377 -0.33026 -0.32636 -0.32208 -0.31741 -0.31234 -0.30686 -0.30097 -0.29466 -0.28793 -0.28076 -0.27315 -0.2651 -0.2566 -0.24764 -0.23822 -0.22832 -0.21794 -0.20708 -0.19573 -0.18387 -0.17151 -0.15863 -0.14523 -0.13131 -0.11684 -0.10184 -0.086285 -0.070173 -0.053496 -0.036246 -0.018417
误差R值：
-0.0089347 -0.01787 -0.026808 -0.035748 -0.044692 -0.05364 -0.062593 -0.071553 -0.08052 -0.089495 -0.098479 -0.10747 -0.11648 -0.12549 -0.13452 -0.14356 -0.15262 -0.16169 -0.17078 -0.17989 -0.18901 -0.19815 -0.20731 -0.21649 -0.2257 -0.23493 -0.24418 -0.25345 -0.26275 -0.27208 -0.28143 -0.29081 -0.30022 -0.30966 -0.31913 -0.32864 -0.33817 -0.34774 -0.35735 -0.36699 -0.37667 -0.38638 -0.39614 -0.40593 -0.41576 -0.42564 -0.43556 -0.44552 -0.45553 -0.46558 -0.47568 -0.48582 -0.49602 -0.50626 -0.51656 -0.5269 -0.5373 -0.54775 -0.55826 -0.56882 -0.57944 -0.59012 -0.60086 -0.61166 -0.62252 -0.63344 -0.64442 -0.65547 -0.66658 -0.67776 -0.68901 -0.70033 -0.71172 -0.72318 -0.73471 -0.74631 -0.75799 -0.76975 -0.78158 -0.79349 -0.80548 -0.81755 -0.8297 -0.84193 -0.85425 -0.86666 -0.87915 -0.89173 -0.9044 -0.91715 -0.93 -0.94295 -0.95599 -0.96912 -0.98235 -0.99568 -1.0091 -1.0226 -1.0363
4. 利用向前差分格式计算一下定解问题

解：根据算法思想，向前差分格式即为
，根据题意，，
，我们可以写出如下程序：
N=10; %把X值从0-1 N等分
M=20; %把Y值从0-1M等分
h=1/N; %X的步长
t=1/M; %Y的步长
X=zeros(1,N+1);
Y=zeros(1,M+1);
Y=0:t:1
X=0:h:1;
r=2*t/(h^2);
U=zeros(M+1,N+1);
%W=zeros(N-1,1);
%S=zeros(N-1,N-1);
%A=zeros(N-1,N-1);
%B=zeros(N-1,N-1);
%C=zeros(N-1,N-1);
%I=zeros(N-1,N-1);
T=zeros(M+1,N+1);
R=zeros(M+1,N+1);
%for j=2:N-2
%S(j,j-1)=1;
%S(j,j+1)=1;
% I(j,j)=1;
%end
%I(1,1)=1;
%I(N-1,N-1)=1;
%S(1,2)=1;
%S(N-1,N-2)=1;
%A=(1-2*r)*I+r*S;
for j=1:N+1 %Y值固定时，X的变化
U(1,j)=exp(X(j))*sin(1/2);
T(1,j)=exp(X(j))*sin(1/2);
end
for j=2:M+1 %X值固定时，Y的变化
U(j,N+1)=exp(1)*sin(1/2-Y(j));
U(j,1)=sin(1/2-Y(j));
T(j,N+1)=exp(1)*sin(1/2-Y(j));
T(j,1)=sin(1/2-Y(j));
end
for k=1:M
for i=2:N
U(k+1,i)=r*U(k,i-1)+(1-2*r)*U(k,i)+r*U(k,i+1); %向前差分格式
T(k+1,i)=exp(X(i))*sin(1/2-Y(k+1));
end
end
for j=1:M+1
for i=1:N+1
R(j,i)=T(j,i)-U(j,i); %误差公式
end
end
运行结果：
真实值T
0.47943
0.52985
0.58557
0.64716
0.71522
0.79044
0.87357
0.96544
1.067
1.1792
1.3032
0.43497
0.48071
0.53127
0.58714
0.64889
0.71714
0.79256
0.87591
0.96803
1.0698
1.1824
0.38942
0.43037
0.47564
0.52566
0.58094
0.64204
0.70957
0.78419
0.86667
0.95781
1.0585
0.3429
0.37896
0.41882
0.46286
0.51154
0.56534
0.6248
0.69051
0.76313
0.84339
0.93209
0.29552
0.3266
0.36095
0.39891
0.44086
0.48723
0.53847
0.5951
0.65769
0.72686
0.80331
0.2474
0.27342
0.30218
0.33396
0.36908
0.4079
0.4508
0.49821
0.55061
0.60852
0.67251
0.19867
0.21956
0.24266
0.26818
0.29638
0.32755
0.362
0.40007
0.44215
0.48865
0.54004
0.14944
0.16515
0.18252
0.20172
0.22294
0.24638
0.27229
0.30093
0.33258
0.36756
0.40621
0.099833
0.11033
0.12194
0.13476
0.14893
0.1646
0.18191
0.20104
0.22218
0.24555
0.27138
0.049979
0.055236
0.061045
0.067465
0.07456
0.082402
0.091068
0.10065
0.11123
0.12293
0.13586
0
0
0
0
0
0
0
0
0
0
0
-0.049979
-0.055236
-0.061045
-0.067465
-0.07456
-0.082402
-0.091068
-0.10065
-0.11123
-0.12293
-0.13586
-0.099833
-0.11033
-0.12194
-0.13476
-0.14893
-0.1646
-0.18191
-0.20104
-0.22218
-0.24555
-0.27138
-0.14944
-0.16515
-0.18252
-0.20172
-0.22294
-0.24638
-0.27229
-0.30093
-0.33258
-0.36756
-0.40621
-0.19867
-0.21956
-0.24266
-0.26818
-0.29638
-0.32755
-0.362
-0.40007
-0.44215
-0.48865
-0.54004
-0.2474
-0.27342
-0.30218
-0.33396
-0.36908
-0.4079
-0.4508
-0.49821
-0.55061
-0.60852
-0.67251
-0.29552
-0.3266
-0.36095
-0.39891
-0.44086
-0.