数值分析常微分数值解的求法C语言

微分方程数值解法C语言-课程设计微分方程数值解法C语言由于对matlab语言不熟悉，所以还是采用C。

前面几个都比较简单，最后一个需要解非其次方程组。

采用高斯—Jordan消元法（数值分析）求逆解方程组，也再一次体会到算法本身的重要性，而不是语言。

当然，矩阵求逆的算法也在100个经典的C语言算法之列。

不过偏微分方程数值解的内容的确比较高深，我只能停留在编这种低级的东西的自娱自乐中。

不过解决计算机、数学、信计专业的课程设计还是足够了。

由于篇幅所限，只把源代码粘贴在这。

一：预报矫正格式#include <math.h>#include<iostream>#include<stdlib.h>double count_0( double xn,double yn){//矫正格式double s;s=yn+0.1*(yn/xn*0.5+xn*xn/yn*0.5);return s;}double count_1(double xn,double yn,double y0){//预报格式double s;s=yn+0.05*((yn/xn*0.5+xn*xn/yn*0.5)+(y0/xn*0.5+xn*xn/y0*0.5));return s;}void main(){//计算,步长为0.1,进行10次计算,设初始值double xn=1,yn=1;int i=1;while(i<=10){printf("%16f ,%1.16f ,%1.16f\n",xn,yn,count_1(xn,yn,count_0(xn,yn)));xn=xn+0.1;yn=count_1(xn,yn,count_0(xn,yn));i++;}}二显示差分格式#include<iostream>#include<math.h>#include<stdlib.h>main(){double a[6][11];//初始化；for(int i=0;i<=5;i++){a[0]=0;a[10]=0;}for(int j=1;j<10;j++){double p=3.14*j*0.1;a[0][j]=sin(p);}//按显示格式计算for(i=1;i<=5;i++)for(j=1;j<10;j++)a[j]=a[i-1][j-1]+a[i-1][j+1]; //输出计算好的矩阵for(i=0;i<=5;i++){for(j=0;j<11;j++)printf("%1.10f ",a[j]);printf("\n");}}三龙阁库塔格式#include <math.h>#include<iostream>#include<stdlib.h>double count_k( double xn,double yn){ double s;s=yn/xn*0.5+xn*xn/yn*0.5;return s;}void main(){//步长为0.1double xn=1,yn=1;int i=1;while(i<=11){printf("%f ,%f\n",xn,yn);double k1=count_k(xn,yn);double k2=count_k(xn+0.05,yn+0.05*k1); double k3=count_k(xn+0.05,yn+0.05*k2); double k4=count_k(xn+0.01,yn+0.1*k3); yn=yn+0.1/6*(k1+2*k2+2*k3+k4);xn=xn+0.1;i++;}}四 CRANK--NICOLSON隐式格式#include<iostream>#include<math.h>#include<stdlib.h>double Surplus(double A[],int m,int n);double * MatrixInver(double A[],int m,int n);double * MatrixOpp(double A[],int m,int n) /*矩阵求逆*/ {int i,j,x,y,k;double *SP=NULL,*AB=NULL,*B=NULL,X,*C;SP=(double *)malloc(m*n*sizeof(double));AB=(double *)malloc(m*n*sizeof(double));B=(double *)malloc(m*n*sizeof(double));X=Surplus(A,m,n);X=1/X;for(i=0;i<m;i++)for(j=0;j<n;j++){for(k=0;k<m*n;k++)B[k]=A[k];{for(x=0;x<n;x++)B[i*n+x]=0;for(y=0;y<m;y++)B[m*y+j]=0;B[i*n+j]=1;SP[i*n+j]=Surplus(B,m,n);AB[i*n+j]=X*SP[i*n+j];}}C=MatrixInver(AB,m,n);return C;}double * MatrixInver(double A[],int m,int n) /*矩阵转置*/ {int i,j;double *B=NULL;B=(double *)malloc(m*n*sizeof(double));for(i=0;i<n;i++)for(j=0;j<m;j++)B[i*m+j]=A[j*n+i];return B;}double Surplus(double A[],int m,int n) /*求矩阵行列式*/ {int i,j,k,p,r;double X,temp=1,temp1=1,s=0,s1=0;if(n==2){for(i=0;i<m;i++)for(j=0;j<n;j++)if((i+j)%2) temp1*=A[i*n+j]; else temp*=A[i*n+j];X=temp-temp1;}else{for(k=0;k<n;k++){for(i=0,j=k;i<m,j<n;i++,j++) temp*=A[i*n+j];if(m-i){for(p=m-i,r=m-1;p>0;p--,r--) temp*=A[r*n+p-1];}s+=temp;temp=1;}for(k=n-1;k>=0;k--){for(i=0,j=k;i<m,j>=0;i++,j--) temp1*=A[i*n+j];if(m-i){for(p=m-1,r=i;r<m;p--,r++) temp1*=A[r*n+p];}s1+=temp1;temp1=1;}X=s-s1;}return X;}void initmat_A(double a[][9],double r){ for(int i=0;i<9;i++)for(int j=0;j<9;j++)a[j]=0;for(i=0;i<9;i++){a=1+r;if(i!=8) a[i+1]=-0.5*r;if(i!=0) a[i-1]=-0.5*r;}}void initmat_B(double b[][9],double r){ for(int i=0;i<9;i++)for(int j=0;j<9;j++)b[j]=0;for( i=0;i<9;i++){b=1-r;if(i!=8) b[i+1]=0.5*r;if(i!=0) b[i-1]=0.5*r;}}void initmat_C(double C[][9]){ for(int i=0;i<9;i++)for(int j=0;j<9;j++)C[j]=0;}void main(){double a[100][11];for(int i=0;i<100;i++)for(int j=0;j<11;j++)a[j]=0;//初始化；for(i=0;i<100;i++){a[0]=0;a[10]=0;}for(int j=1;j<10;j++){double p=4*3.14*j*0.1;a[0][j]=sin(p);}//取h=0.1*3.14,r=0.0005，t=0.0001*3.14*3.14; //得到矩阵a和矩阵bdouble A[9][9];initmat_A(A,0.005);double B[9][9];initmat_B(B,0.005);//B矩阵与Un相乘,en是0;double C[9][9];initmat_C(C);double *A_;A_=MatrixOpp(A[0],9,9);//A矩阵求逆；//A逆*Bfor(i=0;i<9;i++)for(j=0;j<9;j++)for(int s=0;s<9;s++)C[j]+=A_[i*9+s]*B[s][j];//填写a表格for(i=0;i<100;i++){for(j=1;j<10;j++)for(int s=0;s<9;s++)a[i+1][j]+=a[s+1]*C[j-1][s];}//输出表格for(i=0;i<100;i++){for(j=0;j<11;j++)printf("%1.8f ",a[j]);printf("\n");}printf("\n"); printf("\n");//利用精确解，求出表格for(i=0;i<100;i++){for(j=0;j<11;j++)printf("%1.8f",exp(-16*0.0001*0.0005*3.14*3.14*i)*sin(4*j*0.1*3.14));printf("\n");}}。

本科生课程设计报告实习课程数值分析学院名称管理科学学院专业名称信息与计算科学学生姓名学生学号指导教师实验地点实验成绩二〇一六年六月二〇一六年六摘要常微分方程数值解法是计算数学的一个分支.是解常微分方程各类定解问题的数值方法.现有的解析方法只能用于求解一些特殊类型的定解问题，实用上许多很有价值的常微分方程的解不能用初等函数来表示，常常需要求其数值解.所谓数值解，是指在求解区间内一系列离散点处给出真解的近似值.这就促成了数值方法的产生与发展.?关键词：数值解法；常微分1. 实验内容和要求常微分方程初值问题 有精确解2()cos(2)x y x x e x -=+。

要求：分别取步长h = 0.1，0.01，0.001，采用改进的Euler 方法、4阶经典龙格-库塔R －K 方法和4阶Adams 预测-校正方法计算初值问题。

以表格形式列出10个等距节点上的计算值和精确值，并比较他们的计算精确度。

其中多步法需要的初值由经典R-K 法提供。

运算时，取足以表示计算精度的有效位。

2. 算法说明2.1函数及变量说明表1 函数及变量定义1、欧拉方法：1()()(,())i i i i y x y x hf x y x +=+ （i=0，1，2，3，......n-1）(0)y a= (其中a 为初值)2、改进欧拉方法：~1~111()()(,())()()[(,())(,())]2(0)i i i i i i i i i i y x y x hf x y x hy x y x f x y x f x y x y a ++++=+=++=（i=0，1，2......n-1） (其中a 为初值)3、经典K-R 方法： 11213243()6(,)(,)22(,)22(,)i i i i i i i i i i h y y K f x y h hK f x y K h h K f x y K K f x h y hK +⎧=+⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩4、4阶adams 预测-校正方法 预测： 校正：Adsms 内插外插公式联合使用称为Adams 预测-校正系统，利用外插公式计算预测，用内插公式进行校正。

一阶常微分方程数值解的C语言编程实现现实问题与分析学科自身的发展，使得微积分及其涉及领域内出现了寻求数值解的诸多问题(在其他数学领域亦是如此)。

这是由于多种原因造成的，比如现实问题中函数表达式往往并不存在，即使可以拟合也并没有得到表达式的必要且表达式的讨论可能出现较大误差。

再比如，很多复杂函数或函数约束条件很难使用甚至不能使用一般的分析工具，像微分方程存在不可积、不可解的大量实例。

所以数值解的相关理论不断产生、完善。

特别是现代以来，计算机的迅速发展为解决数值解问题提供了算法的计算机程序实现可能，并提供了更多的求解思路。

Mathematica，MATLAB，Maple是最常用的数学软件，对应于Mathematica语言，MATLAB语言，Maple语言。

当然，一般使用C语言也可以比较便捷地解决很多分析学科问题的数值解。

举一元函数微积分学一些熟知的数值解例子：数值导数常见有图解微分法、差商公式、插值多项式、拉格朗日公式(拉格朗日插值公式得来)、马尔科夫公式(牛顿插值公式得来)、等距公式、三次样条函数等等，数值积分有内插求积公式(最为基础的有梯形法、辛普森公式(抛物线法))、高斯型求积公式等等。

这里，考察一阶导数已解出的一阶常微分方程的初值问题的数值解。

问题给出：常微分方程数值解问题中，一阶常微分方程与一阶常微分方程组(一阶导数均解出形式)的数值解法是最基础的，也是非常常用的。

一般数值解法有欧拉方法及改进的欧拉方法、龙格-库塔方法、阿达姆斯方法。

现给出前两种方法，它们都是单步法，而阿达姆斯法是线性多步法，各有优缺点(此外计算机上常会加以使用浮动步长法，如自动变步长计算方法)。

单步法的通用表达式是拉格朗日中值定理的形式(差分代替微分是一种重要的近似思想，反过来用微分代替差分亦重要)：欧拉方法：改进的欧拉方法(预报校正法)：龙格-库塔方法：二阶龙格-库塔方法：四阶龙格-库塔方法：容易看出，欧拉方法及改进的欧拉方法都是特殊的龙格-库塔方法。

使用C语言解常微分方程CODE在C语言中，我们可以使用数值方法来解常微分方程（ODEs）。

常见的数值方法有欧拉法、改进的欧拉法和四阶龙格-库塔法等。

首先，我们需要了解什么是常微分方程。

常微分方程是描述未知函数与其导数之间关系的方程。

一阶常微分方程可以写成如下形式：dy / dx = f(x, y)其中，y是未知函数，f(x,y)表示函数y和自变量x之间的关系。

我们可以通过离散化自变量x的值，来近似求解上述的常微分方程。

假设我们将自变量区间[a,b]划分成N个子区间，每个子区间的长度为h=(b-a)/N，那么我们可以将x离散化为{x0,x1,...,xN}，其中x0=a，xN=b。

对于欧拉法来说，它是最简单的数值方法。

它的基本思想是通过线性逼近来求解常微分方程。

根据导数的定义，可以得到一个线性逼近的公式：dy / dx ≈ (y(i+1) - y(i)) / h将上述式子代入微分方程可以得到：y(i+1)=y(i)+h*f(x(i),y(i))其中x(i+1)=x(i)+h。

下面是一个用C语言实现欧拉法求解一阶常微分方程的例子：```c#include <stdio.h>double f(double x, double y)return x + y;void eulerMethod(double a, double b, int N) double h = (b - a) / N;double x = a;double y = 0; // initial conditionfor (int i = 0; i < N; i++)y=y+h*f(x,y);x=x+h;printf("x = %f, y = %f\n", x, y);}int maidouble a = 0; // initial value of xdouble b = 1; // final value of xint N = 10; // number of intervals eulerMethod(a, b, N);return 0;```在上面的例子中，我们定义了一个常微分方程f(x,y)=x+y。

《计算机数学基础》数值部分第五单元辅导14 常微分方程的数值解法一、重点内容 1. 欧拉公式：),...,,,(),()(1-210=⎩⎨⎧+=+=≈01+1+n k kh x x y x hf y y x y kk k k k k局部截断误差是O (h 2)。

2. 改进欧拉公式：预报－校正公式：⎪⎩⎪⎨⎧++=+=++++)],(),([2),(1111k k k k k k k k k k y x f y x f hy y y x hf y y 校正值预报值即 ))],(,(),([211k k k k k k k k y x hf y x f y x f hy y +++=++ 或表成平均的形式：⎪⎪⎪⎩⎪⎪⎪⎨⎧+21=+=+=1+1+)(),(),(c p k p k k c k k k p y y y y x hf y y y x hf y y改进欧拉法的局部截断误差是O (h 3)3. 龙格－库塔法二阶龙格－库塔法的局部截断误差是O (h 3) 三阶龙格－库塔法的局部截断误差是O (h 4) 四阶龙格−库塔法公式： )22(643211κκκκ++++=+hy y k k其中 κ1=f (x k ,y k )；κ2=f (x n +12h ，y k +21h κ1)；κ3=f (x k +12h ，y n +21h κ2)；κ4=f (x k +h ，y k +h κ3)四阶龙格－库塔法的局部截断误差是O (h 5)。

二、实例例1 用欧拉法解初值问题⎩⎨⎧1=060≤≤0--='2)().(y x xy y y ，取步长h =0.2。

计算过程保留4位小数。

解h =0.2, f (x )=－y －xy 2。

首先建立欧拉迭代格式),,)((.),(210=-420=--=+=21+k y x y y hx hy y y x hf y y k k k kk k k k k k k当k =0，x 1=0.2时，已知x 0=0,y 0=1，有y (0.2)≈y 1=0.2×1(4－0×1)＝0.8000当k ＝1，x 2=0.4时，已知x 1=0.2, y 1=0.8，有 y (0.4)≈y 2=0.2×0.8×(4－0.2×0.8)＝0.614 4 当k =2，x 3=0.6时，已知x 2=0.4,y 2=0.6144，有 y (0.6)≈y 3=0.2×0.6144×(4－0.4×0.4613)＝0.8000例2 用欧拉预报－校正公式求解初值问题⎩⎨⎧1=10=++'2)(sin y x y y y ，取步长h =0.2,计算y (0.2),y (0.4)的近似值，计算过程保留5位小数。

实验七用C++解常微分方程1. 引言常微分方程是自变量和未知函数以及函数的导数之间的关系方程。

在科学领域中，常微分方程被广泛应用于描述物理、生物、经济等领域的各种动力学问题。

在本实验中，我们将使用C++编程语言来解决常微分方程问题。

2. 背景知识在开始解常微分方程之前，我们需要了解几个基本概念和方法：2.1 Euler方法Euler方法是一种基本的数值解常微分方程的方法。

它通过将微分方程转化为差分方程的形式，然后利用差分逼近来递推求解。

其基本公式如下：y_{n+1} = y_n + h*f(x_n, y_n)其中，`y_n` 表示在 `x = x_n` 时的未知函数的值，`f(x_n, y_n)`表示微分方程右端的函数，`h` 表示每一步的步长。

2.2 Runge-Kutta方法Runge-Kutta方法是另一种常用的数值解常微分方程的方法。

相对于Euler方法，Runge-Kutta方法的精度更高。

其中最常用的是4阶Runge-Kutta方法，其基本公式如下：y_{n+1} = y_n + (h/6)*(k_1 + 2k_2 + 2k_3 + k_4)其中，`k_1 = f(x_n, y_n)`，`k_2 = f(x_n + h/2, y_n + (h/2)*k_1)`，`k_3 = f(x_n + h/2, y_n + (h/2)*k_2)`，`k_4 = f(x_n + h, y_n + h*k_3)`。

3. 实验步骤1. 首先，我们需要根据给定的常微分方程，编写一个C++函数用于计算右端函数 `f(x, y)` 的值。

2. 接下来，我们可以选择使用Euler方法或Runge-Kutta方法来求解微分方程。

我们需要编写一个C++函数来实现所选择的方法。

3. 在主函数中，我们可以设定初始条件、步长等参数，并调用函数来求解常微分方程。

4. 最后，将求得的数值解绘制成图表，以便对结果进行可视化分析。

i.常微分方程初值问题数值解法i.1 常微分方程差分法考虑常微分方程初值问题：求函数()u t 满足(,), 0du f t u t T dt=<≤ （i.1a ） 0(0)u u = (i.1b)其中(,)f t u 是定义在区域G : 0t T ≤≤, u <∞上的函数，0u 和T 是给定的常数。

我们假设(,)f t u 对u 满足Lipschitz 条件，即存在常数L 使得121212(,)(,), [0,]; ,(,)f t u f t u L u u t T u u -≤-∀∈∈-∞∞ （i.2） 这一条件保证了（i.1）的解是适定的，即存在，唯一，而且连续依赖于初值0u 。

通常情况下，（i.1）的精确解不可能用简单的解析表达式给出，只能求近似解。

本章讨论常微分方程最常用的近似数值解法--差分方法。

先来讨论最简单的Euler 法。

为此，首先将求解区域[0,]T 离散化为若干个离散点：0110N N t t t t T -=<<<<= （i.3） 其中n t hn =，0h >称为步长。

在微积分课程中我们熟知，微商（即导数）是差商的极限。

反过来，差商就是微商的近似。

在0t t =处，在（i.1a ）中用向前差商10()()u t u t h -代替微商du dt ，便得 10000()()(,())u t u t hf t u t ε=++如果忽略误差项0ε，再换个记号，用i u 代替()i u t 便得到1000(,)u u hf t u -=一般地，我们有1Euler (,), 0,1,,1n n n n u u hf t u n N +=+=-方法： （i.4） 从(i.1b) 给出的初始值0u 出发，由上式可以依次算出1,,N t t 上的差分解1,,N u u 。

下面我们用数值积分法重新导出 Euler 法以及其它几种方法。

为此，在区间1[,]n n t t +上积分常微分方程（i.1a ），得11()()(,())n n t n n t u t u t f t u t dt ++=+⎰ （i.5）用各种数值积分公式计算（i.5）中的积分，便导致各种不同的差分法。