?

#include <stdio.h>

#include <cmath>

#define N 20//離散點個數取20，對應步長0.1

#define low -1.0//定義域x區間下限

#define high? 1.0//定義域x區間上限

#define h (high-low)/N//步長

double q(double x)//定義函數q(x)

{

????? return 1+pow(x,2);

}

double gamma(double x)//定義函數gamma(x)

{

????? return 0;

}

int main()

{

????? double x[N+1];//自變量

????? double y[N+1];//y

????? double alpha=1,beta=1;//邊界值

????? double Gamma[N];

????? double Q[N];

????? double a[N-2];//三對角陣中主對角線下方元素

????? double b[N-1];//三對角陣中主對角線的元素

????? double c[N-2];//三對角陣中主對角線上方元素

????? double f[N-1];

????? double z[N-1];

????? double M[N-1];

????? double Alpha[N-1];

????? double Beta[N-2];

????? x[0]=low, x[N-1]=high;

????? y[0]=1, y[N]=1;

????? int i;//循環計數變量

?????

????? for(i=1;i<N+1;i++)//將區間進行N等分，記分點為x[i]

????? {

???????????? x[i]=x[0]+i*h;

????? }

?????

????? for(i=0;i<N;i++)

????? {

???????????? Q[i]=q(x[i+1]);

????? }

?????

????? for(i=0;i<N;i++)

????? {

???????????? Gamma[i]=gamma(x[i+1]);

????? }

????? //確定近似值y[i]的線性方程組，改寫成系數矩陣為對角占優的矩陣形式（三對角矩陣）

????? //用追趕法解三對角線方程組AM=f,其中M為要解的未知量

????? //給三對角陣A賦值

????? for(i=0;i<N;i++)

????? {

???????????? b[i]=-2-Q[i]*pow(h,2);

????? }

????? for(i=0;i<N-1;i++)

????? {

???????????? a[i]=1;

????? }

????? for(i=0;i<N-1;i++)

????? {

???????????? c[i]=1;

????? }

?????

????? f[0]=Gamma[0]*pow(h,2)-alpha;

????? f[N-2]=Gamma[N-2]*pow(h,2)-beta;

????? for(i=1;i<N-2;i++)

????? {

???????????? f[i]=Gamma[i]*pow(h,2);

????? }

?????

????? //A=LU分解

????? Beta[0]=c[0]/b[0];

????? for(i=1;i<N-2;i++)

????? {

???????????? Beta[i]=c[i]/(b[i]-a[i-1]*Beta[i-1]);

????? }

????? Alpha[0]=b[0];

????? for(i=1;i<N-1;i++)

????? {

???????????? Alpha[i]=b[i]-a[i-1]*c[i-1]/Alpha[i-1];

????? }

????? //解Lz=f

????? z[0]=f[0]/b[0];

????? for(i=1;i<N-1;i++)

????? {

???????????? z[i] = (f[i] - a[i - 1] * z[i - 1]) / Alpha[i];

????? }

????? //解UM=z

????? M[N-2]=z[N-2];

????? for(i=N-3;i>=0;i--)

????? {

???????????? M[i]=z[i]-(Beta[i]*M[i+1]);

????? }

????? for(i=1;i<N;i++)

????? {

???????????? y[i]=M[i-1];

????? }

????? for(i=1;i<N;i++)

????? {

???????????? printf("y(%0.1f)=%0.6f\n",x[i],y[i]);

????? }

????? return 0;

}

?

總結

以上是生活随笔為你收集整理的常微分方程数值解法（2）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。