C語言編程解線性,非線性方程,龍貝格算法

1、課 程 設 計課程名稱： 數(shù)值計算B 設計題目： 數(shù)值計算B課程設計 學 號： 姓 名： 完成時間： 2015年10月27日 題目一：非線性方程求根一、題目假設人口隨時間和當時人口數(shù)目成比例連續(xù)增長，在此假設下人口在短期內的增長建立數(shù)學模型。（1）如果令表示在時刻的人口數(shù)目，表示固定的人口出生率，則人口數(shù)目滿足微分方程，此方程的解為；（2）如果允許移民移入且速率為恒定的，則微分方程變成，此方程的解為；假設某地區(qū)初始有1000000人，在第一年有435000人移入，又假設在第一年年底該地區(qū)人口數(shù)量1564000人，試通過下面的方程確定人口出生率，精確到；且通過這個數(shù)值來預測第二年年末的人口數(shù)，假

2、設移民速度保持不變。二、數(shù)學原理本題考查非線性方程的解法。數(shù)值計算中常用的有二分法、不動點迭代法、牛頓法以及其變形，本文采用牛頓法進行迭代計算。牛頓法的基本思想是將非線性方程逐步轉化為某種線性方程來求解。其計算過程大致為：（1） 將非線性方程轉化為的形式；（2） 將函數(shù)在方程的近似根（假定）處泰勒展開，則原方程可近似為（3） 記上述方程的根為，得到其計算公式為（4）按照給定的精度要求判斷是否收斂，得到最終結果。本題中，將非線性方程轉化為的形式如下：即則由于牛頓法對初值的要求較高，因此先采用二分法初步給定初值。如下表所示：0.11的符號-+ 故取初值。三、程序設計使用C語言進行迭代程序的設計。定

3、義數(shù)組x1000存放人口增長率的迭代數(shù)據(jù)；定義變量a表示迭代次數(shù)，定義變量b為迭代誤差，變量c為計算精度；輸入變量：x0=0.5為初值，c=0.0001為計算精度；輸出變量：xa（即最終結果），a，迭代總次數(shù)，b迭代最終誤差。程序設計見附錄。（二分法與牛頓法混合使用，計算收斂較易，因而只給定最大迭代次數(shù)，未進行迭代是否發(fā)散的判定）四、結果分析和討論程序運行結果如下圖所示，計算精度為10-4的最終結果為。將其代入方程得到第二年人口的計算式為即人口出生率為10.10%，第二年該地區(qū)人口數(shù)量約為2187900人。討論：使用牛頓法與二分法相配合在較少的迭代次數(shù)下即取得了較為準確的結果，可見數(shù)值計算方法

4、在計算類似問題上有著很大的優(yōu)勢。圖1 題目一程序運行結果五、完成題目的體會與收獲在求解該題目的過程中，使用了非線性方程的牛頓迭代法求解，運用C語言設計了迭代計算程序，對非線性方程求解的數(shù)學思想有了深入認識和切實體會，有助于將數(shù)值計算課程中學到的理論知識運用到實際物理問題的研究中。題目二：線性方程組求解一、題目假設一個物體可以位于個等距點的任意位置，當物體在位置時，它只能等可能的移動到或者，而不能直接移動到其他任何位置，概率表示物體從位置開始在到達右端點之前到達左端點的概率，顯然，且有，即有下面方程組：取對方程組進行求解（迭代法或者直接法）。二、數(shù)學原理觀察上述線性方程組的系數(shù)矩陣，顯然該矩陣是

5、弱對角占優(yōu)的三對角陣，可利用三對角方程組的追趕法求解。給出三對角陣（針對本題系數(shù)矩陣，給出九階方陣克勞特分解）的克勞特分解如下： 其中為待定系數(shù)，追趕法公式及計算過程如下：（1） 計算待定系數(shù)的遞推公式： 得（2） 求解，（3）求解： 由此設計計算程序，并輸入系數(shù)即可完成線性方程組的求解。三、程序設計使用C語言進行計算程序的設計。定義數(shù)組a8,b9,c8存放三對角系數(shù)矩陣的三類值；定義數(shù)組d9存放方程組右側向量；定義m8存放待定系數(shù)的值；定義數(shù)組x9,y9存放方程組的解。輸入變量：a8,b9,c8,d9；輸出變量：x9。程序設計見附錄。四、結果分析和討論程序運行結果如下圖所示：圖2 題目二運算

6、結果即原方程組的解為考慮，并將方程組的解中的無效數(shù)位去掉，則原物理問題的解向量（11階）為 討論：使用追趕法直接求解三對角占優(yōu)的系數(shù)矩陣所在的線性方程組計算結果與迭代方法相比在原理上更為可靠，且充分利用了系數(shù)矩陣的特點，大大減小了計算量和存儲空間，在解決類似問題的過程中優(yōu)勢明顯。五、完成題目的體會與收獲通過本題目的計算，深入體會了追趕法求解三對角系數(shù)矩陣線性方程組的“追趕”過程，在實踐中加強了對理論知識的理解。 題目三：數(shù)值積分一、題目衛(wèi)星軌道是一個橢圓，橢圓周長的計算公式是, 這里是橢圓的半長軸, 是地球中心與軌道中心(橢圓中心)的距離, 記為近地點距離, 為遠地點距離, 公里為地球半徑,則

7、, 某人造衛(wèi)星近地點距離公里,遠地點距離公里, 試用Romberg方法求衛(wèi)星軌道的周長，精確到。二、數(shù)學原理龍貝格方法在思想上不同于插值型求積公式，其目的是對近似值進行修正而得到更精確的計算公式，實質上是對復化梯形公式的多次修正。以下簡要介紹龍貝格算法的計算過程。待求解定積分為，以表示二分k次后求得的梯形值，以表示序列的m次加速值，理查森外推法的遞推公式可寫成計算過程如下：（1） 取，求；（2） 利用變步長梯形公式求，其中k為區(qū)間的二分次數(shù)，即（3）橫向求解；（4）相鄰兩次計算的小于給定的精度要求時停止計算。本題中根據(jù)已知條件，積分上下限分別為、0。積分中涉及到的常數(shù)計算如下：被積函數(shù)為三、程

8、序設計使用C語言進行計算程序的設計。定義變量a，b存放積分下限與上限；定義數(shù)組T100100存放迭代結果，Tii即第i次迭代的積分值；定義變量n存放誤差限；定義變量m，w存放最終迭代誤差和總迭代次數(shù)。輸入變量：；輸出變量：計算結果Tww，迭代次數(shù)w，最終誤差m。 程序設計見附錄。四、結果分析和討論程序計算結果如下圖所示：圖3 問題三程序結果即公里公里 1/4軌道周長為12331.298828公里，整個軌道周長為49325.195312公里。討論：（1）由于積分上限中含有，其值難以確定，在程序中提前賦值，保留10位小數(shù)，由此引入部分誤差；（2） 題中給出的計算公式應為橢圓1/4軌道周長，整個軌道

9、周長應為4a。且由于計算S產(chǎn)生的誤差遠小于10-6，則4S的誤差也小于10-6，即4S的計算精度也在題目要求的范圍內。五、完成題目的體會與收獲使用龍貝格算法完成了數(shù)值積分的計算，迭代4次即得到了精度為10-6的計算結果，充分展示了數(shù)值計算的重要作用，同時有助于領會龍貝格算法中體現(xiàn)的數(shù)學思想，為解決今后可能遇到的物理問題提供了保證。附錄題目一#include<stdio.h>#include<stdlib.h>#include<string.h>#include<math.h>void main()/定義變量double x1000;/存放人口出生

10、率迭代結果，限定最大迭代次數(shù)為1000次x0=0.5;/為人口出生率賦初值double f,f1;/存放函數(shù)的原函數(shù)及導函數(shù)迭代結果int a;/迭代次數(shù)double b;/相鄰兩次迭代間誤差double c=0.0001;/誤差限/求解過程for(int i=0;i<1000;) f=1000000*exp(xi)+435000*(exp(xi)-1)/xi-1564000; f1=1000000*exp(xi)+(435000*exp(xi)/xi-(435000*exp(xi)-435000)/pow(xi,2); xi+1=xi-f/f1; b=fabs(xi+1-xi); pr

11、intf("第%d次迭代計算誤差: %.5fn",i+1,b); if(b<c) a=i+1;break; else i+; /結果輸出printf("n計算結果: %.4f",xa);printf("n總迭代次數(shù): %d",a);printf("n計算誤差: %.6f",b);system("pause");題目二#include<stdio.h>#include<stdlib.h>#include<string.h>#include<math.

12、h>void main()/定義變量double a8=-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5;double b9=1,1,1,1,1,1,1,1,1;double c8=-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5,-0.5;/存放系數(shù)三對角矩陣系數(shù)double d9=0.5,0,0,0,0,0,0,0,0;/存放方程組右側數(shù)據(jù)double m8;/存放待定系數(shù)貝塔double x9,y9;/存放方程組的解/求解過程m0=c0/b0;for(int i=1;i<8;i+)mi=ci/(bi-ai-1*mi-1);y0

13、=d0/b0;for(int j=1;j<9;j+)yj=(dj-aj-1*yj-1)/(bj-aj-1*mj-1);x8=y8;for(int k=7;k>=0;k-)xk=yk-mk*xk+1;/結果輸出printf("計算結果為n");for(int i1=0;i1<9;i1+)printf("x%d=%.4fn",i1+1,xi1);system("pause");附錄三#include<stdio.h>#include<stdlib.h>#include<string.h>

14、;#include<math.h>float f(float x)return(sqrt(7880.5*7880.5-973.5*973.5*sin(x)*sin(x);/被積函數(shù)void main()/定義變量float pi=3.1415926535;/為賦值float a=0,b=0.5*pi;/定義積分上下限float T100100;/存放迭代結果T00=(b-a)/2*(f(a)+f(b);float m;/存放歷次迭代誤差float n=0.000001;/誤差限int w;/迭代次數(shù)float ys;float y;/計算過程for(int i=1;i<100;i+) int p=pow(float)2,i-1)-1; ys=0;for(int k=0;k<=p;k+) y=f(a+(2*k+1)*(b-a)/pow(float)2,i); ys=ys+y; Ti0=0.5*Ti-10+(b-a)/pow(float)2,i)*ys;for(int j=1;j<=i;j+)Tij=pow(float)4,j)/