平淡才真i
1、 编程实现以下科学计算算法,并举一例应用之。 “中点公式法和五点公式法求数值微分”
解:中点公式法和五点公式法求数值微分如下:
例 5-4:中点公式法求导数应用实例。 采用中点公式法求函数 f=x 在 x=4处的导数。
解:在 MATLAB 命令窗口中输入:
>>df=MidPoint('sqrt(x)',4)
输出结果为:
df=
0.2500
采用中点公式法求函数 f=x 在 x=4处的导数为 0.25,而导数的精确值也是 0.25 .
详见以下:
中点公式法流程图:
源代码:
function df=MidPoint(func,x0,h)
if nargin == 2
h = 0.1;
else if (nargin == 3 && h == 0.0)
disp('h2??ü?a0￡?');
return ;
end
end
y1 = subs(sym(func), findsym(sym(func)),x0+h);
y2 = subs(sym(func), findsym(sym(func)),x0-h);
df = (y1-y2)/(2*h);
运行结果如下:
例 5-5:五点公式法求导数应用实例。 采用五点公式法求函数 f=sin(x)在 x=2处的导数。
解:在 MATLAB 命令窗口中输入:
>>df1=FivePoint('sin(x)',2,1);
>>df2=FivePoint('sin(x)',2,2);
>>df3=FivePoint('sin(x)',2,3);
>>df4=FivePoint('sin(x)',2,4);
>>df5=FivePoint('sin(x)',2,5);
用五种方法得到的结果为:
df1=
-0.4161
df2=
-0.4161
df3=
-0.4161
df4=
-0.4161
df5=
-0.4161
而函数在 f=sin(x)在 x=2的导数为 cos(2)=-0.4161,从上面的结果来看,五点公式的 精度是很高的。
详见以下:
五点公式法流程图 :
源代码:
function df=FivePoint(func,x0,type,h)
%??μ?1?ê?·¨,?óè?oˉêyfunc?úx0′|μ?êy ×￡íòê?è?òa
%oˉêy??￡ofunc
%?óμ?μ?￡ox0
%1?ê?μ?D?ê?￡otype￡¨è?1,2,3,4,5,￡?
%à?é￠2?3¤￡oh
%μ?êy?μ￡odf
if nargin ==3
h=0.1;
else if (nargin ==4&&h==0.0)
disp('h2??ü?a0');
return ;
end
end
y0 = subs(sym(func),findsym(sym(func)),x0);
y1 = subs(sym(func),findsym(sym(func)),x0+h);
y2 = subs(sym(func),findsym(sym(func)),x0+2*h);
y3 = subs(sym(func),findsym(sym(func)),x0+3*h);
y4 = subs(sym(func),findsym(sym(func)),x0+4*h);
y_1 = subs(sym(func),findsym(sym(func)),x0-h);
y_2 = subs(sym(func),findsym(sym(func)),x0-2*h);
y_3 = subs(sym(func),findsym(sym(func)),x0-3*h);
y_4 = subs(sym(func),findsym(sym(func)),x0-4*h);
switch type
case 1,
df=(-25*y0+48*y1-36*y2+16*y3-3*y4)/(12*h);%ó?μúò???1?ê??óμ?êycase 2,
df=(-3*y_1-10*y0+18*y1-6*y2+y3)/(12*h);%ó?μú?t??1?ê??óμ?êy case 3,
df=(y_2-8*y_1+8*y1-y2)/(12*h);%ó?μúèy??1?ê??óμ?êy case 4,
df=(3*y1+10*y0-18*y_1+6*y_2-y_3)/(12*h);%ó?μú????1?ê??óμ?êy case 5,
df=(25*y0-48*y_1+36*y_2-16*y_3+3*y_4)/(12*h);%ó?μú????1?ê??óμ?êy end
运行结果如下:
2、编程解决以下科学计算和工程实际问题。
① 已知阿波罗(Apollo )卫星的运动轨迹 (x,y)满足下列微分方程
()r
r x x x y
x 32
*31
*..
)
(2μμμμ--
+-+=
r
r
y
y
y x y 32
31
*.
..
2μμ-
-+-=
其中 μ=45
. , *
μ=1-μ
221
) (y x r
++=μ , 22*2) (y x r ++=μ 试在初值
x(0)=1.2, 0) 0(.
=x , , 04935751. 1) 0(.
-=y 下进行数值求解, 并绘制出阿波罗卫星位置 (x,y)
的轨迹。
① 解:根据题目选用 MATLAB 代码如下:
function dy=weifen(t,y)
% 编程解决阿波罗(Apollo )卫星的运动轨迹 求解器属于变步长的一种,
采用 Runge-Kutta 算法 万事如意
u=1/82.45; b=1-u;
dy=zeros(4,1); r=zeros(2,1);
r(1)=sqrt((y(1)+u)^2+(y(3))^2); r(2)=sqrt((y(1)+b)^2+(y(3))^2); dy(1)=y(2);
dy(2)=2*dy(3)+y(1)-b*(y(1)+u)/(r(1)^3)-u*(y(1)-b)/(r(2)^3); dy(3)=y(4);
dy(4)=-2*dy(1)+y(3)-b*y(3)/(r(1)^3)-u*y(3)/(r(2)^3);
解:在 MATLAB 命令窗口中输入
>>ode45('weifen',[0 2.00],[1.2 0 0 -1.04935751])
>>[T,Y]=ode45('weifen',[0 1.26],[1.2 0 0 -1.04935751])
运行结果:
阿波罗卫星位置 (x,y)的轨迹图如下:
② 实验图所示是一个跷跷板,两板价较为,左边板长为 1.5m ,上面的小孩重 150N, 右 边板长为 2m, 小孩重为 400N. 求当跷跷板平衡时, 左边木板与水平方向的夹角ɑ的大小。 要求 先求解析解,然后给出两种解决方案。
② 解:根据力矩平衡求解析解
由图示可有下列关系式:
500?1.5αcos =2?400) 3
1cos(απ- 解该式得:
7
87
8tan sin cos 350==
=ααα
α 即:rad 4678. 0≈α
两种方法的求解:
方案一:采用两步迭代法求解方程:
500?1.5αcos =2?400) 3
1
cos(απ-
两步迭代法的 MATLAB 的代码如下:
源代码:
function root=TwoStep(f,a,b,type,eps) if (nargin==4)
eps=1.0e-4; %eps ±íê?μü′ú???è end
f1=subs(sym(f),findsym(sym(f)),a);
f2=subs(sym(f),findsym(sym(f)),b);
if (f1==0)
root=a;
end
if (f2==0)
root=b;
end
if (f1*f2>0)
disp(' 两端点函数值乘积大于 0!' ); return ;
else
tol=1;
fun=diff(sym(f));
fa=subs(sym(f),findsym(sym(f)),a);
fb=subs(sym(f),findsym(sym(f)),b);
dfa=subs(sym(fun),findsym(sym(fun)),a);
dfb=subs(sym(fun),findsym(sym(fun)),b);
if (dfa>dfb)
root=a;
else
root=b;
end
while (tol>eps)
if (type==1)
r1=root;
fx1=subs(sym(f),findsym(sym(f)),r1);
dfx=subs(sym(fun),findsym(sym(fun)),r1);
r2=r1-fx1/dfx;
fx2=subs(sym(f),findsym(sym(f)),r2);
root=r2-fx2/dfx;
tol=abs(root-r1);
else
r1=root;
fx1=subs(sym(f),findsym(sym(f)),r1);
dfx=subs(sym(fun),findsym(sym(fun)),r1);
r2=r1-fx1/dfx;
fx2=subs(sym(f),findsym(sym(f)),r2);
root=r2-fx2*(r2-r1)/(2*fx2-fx1);
tol=abs(root-r1);
end
end
end
解:在 MATLAB 命令窗口中输入
r=TwoStep('500*1.5*cos(x)-2*400*cos(1/3*pi-x)',0,1/3*pi,1) 运行结果如下:
两个小孩所产生力矩随 α变化的曲线如下图片:
运用 Data cursor 工具,得到交点处对应的 X 值为 :0.4678 也即:4678. 0=α
>> x=0:0.0001:1/3*pi;
>> ML=500*1.5*cos(x);
>> MR=400*2*cos(1/3*pi-x); >> plot(x,ML ,'-',x,MR ,'-')
数值分析中的各种公式C 代码
二分法
2.2 算法步骤
步骤 1:准备 计算 f(x)在有根区间 [a,b]端点处的值 f(a), f(b).
步骤 2:二分 计算 f(x)在区间中点 (a+b)/2处的值 f((a+b)/2).
步骤 3:判断 若 f((a+b)/2)=0,则 (a+b)/2即是根,计算过程结束,否则检验;
若 f((a+b)/2)f(a)<0,则以 (a+b)/2代替="" b="" ,否则以="" (a+b)/2代替="">
反复执行步骤 2和步骤 3,直到区间 [a,b]的长度小于允许误差 e ,此时中点 (a+b)/2即为所 求近似根。
2.3 程序流程图
3 实验结果分析
#include
void main()
{
float a,b,e,x,c;
int k=0,n=1;
scanf(
while(n==1)
{
x=(a+b)/2;c=(x*x*x-x-1)*(a*a*a-a-1); if(c<>
{
b=x;
if(b-a<>
{
printf(
else k++;}
else
{ if(c==0) { printf(
break;}
else { a=x; if(b-a<=e) {="">
break;}
else k++;
}}}}
高斯塞德尔迭代法求方程组解 高斯主元素消去法求方程组解
2.2 算法步骤
高斯塞德尔迭代法:
步骤 1:给定初值 )
0(1) 0(2) 0(1,..., , n x x x ,精度 ε,最大容许迭代次数 M, 令 k=1。
步骤 2:对 i=1,2,… ,n 依此计算 )
0() 1()
0()
1(01) 0() 1()
1(). (i i i
i i ii
n
j j j ij i
i
x x x x e a x a x
x
→-=-=∑≠=
步骤 3:求出 e=
}{max 1i n
i e ≤≤,若 e<ε,则输出>
0(i x (i=1,2,..,n ),停止计算。否则执行步 步骤 4:若 k
高斯主元素消去法:
步骤 1:det ←1, k=1,2,… ,n-1, 从第 2步做到第 7步; 步骤 2: 按列主元素 ik k ik a a n
i k , max ≤≤=
步骤 3:如果 k ik a , =0,则 det ←0,计算停止。 步骤 4:若 k i k =,转步骤 5,否则换行:
det det , ), ,... 1, (, , ←?+=?ik k j ik kj b b n k k j a a
步骤 5:计算乘数 ik m , kk ik ik ik a a m a /=← (i=k_1,… ,n) 步骤 6:消元计算
)
,..., 1(, ) ,..., 1, (, n k i b m b b n k j i a m a a i ik i i kj ik ij ij +=-←+=-←
步骤 7:det kk a ←det;
步骤 8:若 0det , 0←=nn a 计算停止。否则 det det nn a ←
步骤 9:回代求解
ii
j
n
i j ij
ij i nn
n n a
b a a b a b b ) (/1
∑+=-
←← (i=n-1,… ,1)
2.3 程序流程图
3 实验结果分析 高斯塞德尔迭代法:
高斯主元素消去法:
高斯塞德尔迭代法:
#include #include #include #define M 10 int n; void Gauss_seidel(double x[],double x0[],double a[][M],double b[],double eps,int Nmax) { int i,j,s=0; double max; while(s<> { for(i=0;i<> { x[i]=b[i]; for(j=0;j<> if(j!=i) x[i]=x[i]-a[i][j]*x[j]; x[i]=x[i]/a[i][i]; } max=fabs(x[0]-x0[0]); for(i=1;i<> if(fabs(x[i]-x0[i])>max) max=fabs(x[i]-x0[i]); if(max<> break; for(i=0;i<> x0[i]=x[i]; s++; } if(s>=Nmax) cout { cout printf( } } void main() { double a[10][10]={{5,2,1},{-1,4,2},{2,-3,10}}; double b[10]={-12,20,3}; double x[10],x01[10]; double eps;// 精度 int i,j,s=0,Nmax; cout n=3; for(i=0;i<> { for(j=0;j<> cout<> cout<><><> } cout for(i=0;i<> { printf( cin>>x[i]; x01[i]=x[i]; } cout cin>>eps; cout cin>>Nmax; Gauss_seidel(x,x01,a,b,eps,Nmax); cout<> } 高斯主元素消去法: #include #include #include //在列向量中寻找绝对值最大的项,并返回该项的标号 int FindMax(int p,int N,double *A) { int i=0,j=0; double max=0.0; for(i=p;i<> { if(fabs(A[i*(N+1)+p])>max) { j=i; max=fabs(A[i*(N+1)+p]); } } return j; } //交换矩阵中的两行 void ExchangeRow(int p,int j,double *A,int N) { int i=0; double C=0.0; for(i=0;i<> { C=A[p*(N+1)+i]; A[p*(N+1)+i]=A[j*(N+1)+i]; A[j*(N+1)+i]=C; } } //上三角变换, A 为增广矩阵的指针, N 为矩阵的行数。 void uptrbk(double *A,int N) { int p=0,k=0,q=0,j=0; double m=0.0; for(p=0;p<> { //找出该列最大项的标号 j=FindMax(p,N,A); //交换 p 行和 j 行 ExchangeRow(p,j,A,N); if(A[p*(N+1)+p]==0) { printf( } //消去 P 元素以下的 p 列内容。 for(k=p+1;k<> { m=A[k*(N+1)+p]/A[p*(N+1)+p]; for(q=p;q<> A[k*(N+1)+q]=A[k*(N+1)+q]-m*A[p*(N+1)+q]; } } printf( { if(j%(N+1)==0) printf( printf( } } //下面是回代函数 double* backsub(double *A,int N) { double* X=NULL,temp=0.0; int k=0,i=0; X=(double*)malloc(N*sizeof(double)); X[N-1]=A[(N-1)*(N+1)+N]/A[(N-1)*(N+1)+N-1]; for(k=N-2;k>=0;k--) { temp=0.0; for(i=k+1;i<> temp=temp+A[k*(N+1)+i]*X[i]; X[k]=(A[k*(N+1)+N]-temp)/A[k*(N+1)+k]; } return X; } main() { int N=0,i=0; double *A=NULL,*X=NULL; printf( scanf( if(N>256||N<> { printf( printf( return 0;} else { A=(double*)calloc(N*(N+1),sizeof(double)); printf( scanf( system( printf( for(i=0;i<> { if(i%(N+1)==0) printf( printf( uptrbk(A,N); //上三角变换 X=backsub(A,N); //回代函数 printf( for(i=0;i<> printf( free(A); free(X); exit(0); } 二次或者三次样条插值 #include #include #include #include #include #define BOOL int #define FALSE 0 #define TRUE 1 //计算分段线性插值的系数 void LinearCoe(double * pd_X, double * pd_Y,double * pd_Coe, int i) { pd_Coe[0] = (-pd_X[i]*pd_Y[i+1]+pd_X[i+1]*pd_Y[i])/(pd_X[i+1]-pd_X[i]); pd_Coe[1] = (pd_Y[i+1]-pd_Y[i])/(pd_X[i+1]-pd_X[i]); } //计算分段二次多项式插值的系数 void QuadrateCoe(double * pd_X, double * pd_Y,double * pd_Coe, int i) { int n = 2; double * pd_MatrixA = new double[(n+1)*(n+1)]; double * pd_VectorB = new double[n+1]; pd_VectorB[0] = pd_Y[i-1]; pd_VectorB[1] = pd_Y[i]; pd_VectorB[2] = pd_Y[i+1]; for(int j = 0; j < n+1;=""> { pd_MatrixA[j*3] = pd_X[i+j-1]*pd_X[i+j-1]; pd_MatrixA[j*3 + 1] = pd_X[i+j-1]; pd_MatrixA[j*3 + 2] = 1; } GaussRemove(pd_MatrixA,pd_VectorB,pd_Coe,n+1); } //计算分段三次多项式插值的系数 void ThriceCoe(double * pd_X, double * pd_Y,double * pd_Coe, int i) { int n = 3; double * pd_MatrixA = new double[(n+1)*(n+1)]; double * pd_VectorB = new double[n+1]; pd_VectorB[0] = pd_Y[i-2]; pd_VectorB[1] = pd_Y[i-1]; pd_VectorB[2] = pd_Y[i]; pd_VectorB[3] = pd_Y[i+1]; for(int j = 0; j < n+1;=""> { pd_MatrixA[j*4] = pd_X[i+j-2]*pd_X[i+j-2]*pd_X[i+j-2]; pd_MatrixA[j*4 + 1] = pd_X[i+j-2]*pd_X[i+j-2]; pd_MatrixA[j*4 + 2] = pd_X[i+j-2]; pd_MatrixA[j*4 + 3] = 1; } GaussRemove(pd_MatrixA,pd_VectorB,pd_Coe,n+1); } //计算三次样条插值函数的系数组 M void SplineCoe(double * pd_X, double * pd_Y, int N, double * pd_h, double * pd_M) { int n = N -1; double * pd_Alpha = new double[n]; //系数数组 Alpha double * pd_Beta = new double[n]; //系数数组 Beta double * pd_Gama = new double[n]; //系数数组 Gama int i,j; for(i = 0; i < n;="" i++)="" 计算=""> { pd_Alpha[i] = pd_h[0]/(pd_h[0] + pd_h[i]); pd_Gama[i] = 1 - pd_Alpha[i]; pd_Beta[i] = 6 * (((pd_Y[1]-pd_Y[0])/pd_h[0]) - ((pd_Y[i+1]-pd_Y[i])/pd_h[i])) / (pd_h[0] + pd_h[i]); } double * pd_MatrixA = new double[n*n]; for(i = 0; i < n;=""> { for(j = 0; j < n;=""> { if (i == j) { pd_MatrixA[i*n+j] = 2; } else if (j == (i+1)) { pd_MatrixA[i*n+j] = pd_Alpha[i]; } else if (j == (i-1)) { pd_MatrixA[i*n+j] = pd_Gama[i]; } else pd_MatrixA[i*n+j] = 0; } } pd_MatrixA[n-1] = pd_Gama[0]; pd_MatrixA[(n-1)*n] = pd_Alpha[n-1]; // ShowMatrix(pd_MatrixA,n,n); GaussRemove(pd_MatrixA,pd_Beta,&(pd_M[1]),N-1); pd_M[0] = pd_M[N-1]; // ShowVector(pd_Beta,n); delete pd_Alpha; delete pd_Beta; delete pd_Gama; delete pd_MatrixA; } 复化梯形公式和复化辛普森公式 对应程序 实验一的源程序: #include #include double fun1(double x,double y,int n) { int i; double sum=0,m; for(i=1;i<=n-1;i++) 求="" f(x(k))的和=""> m=x+i*y; sum=sum+sqrt(1+pow(2.718,m)); } return sum; } double get(double x) { return sqrt(1+pow(2.718,x)); } double fun2(double x,double y,int n) { int i; for(i=1;i<=n-1;i++) 求=""> { sum1+=get(x+i*y+y/2); sum2+=get(x+i*y); } ff=4*sum1+2*sum2; return ff; } void main() { double a=0,b=2;//给定的区间 int n; while(1) { cout double h; double I1,f10,f11,f12; double I2,f20,f21,f22; h=(b-a)/n;//步长 f10=get(a); f11=fun1(a,h,n); f12=get(b); f20=f10; f21=fun2(a,h,n); f22=f12; I1=h/2*(f10+2*f11+f12);//复化梯形公式 I2=h/6*(f20+f21+f22); //复化辛浦生公式 cout cout } 实验内容二的源程序: #include #include double get(double x) { if(x==0) return 1; else return sin(x)/x; } double fun(double x,double y,int n) { int i; for(i=1;i<=n-1;i++) 求=""> { sum1+=get(x+i*y+y/2); sum2+=get(x+i*y); } ff=4*sum1+2*sum2; return ff; } void main() { double a=0,b=1;//给定的区间 int n;//等分成 n 个空间 while(true) { cout cin>>n; if(n==-1) break; double h; double I1,f11,f12,f13; h=(b-a)/n;//步长 f11=get(a); f12=fun(a,h,n); f13=get(b); I1=h/6*(f11+f12+f13);//复化辛浦生公式 cout cout } } 10个重要的算法 C 语言实现源代码:拉格朗日,牛顿插值,高斯,龙贝格, 牛顿迭代,牛顿 -科特斯,雅克比,秦九昭,幂法,高斯塞德尔 (转 ) 1. 拉格朗日插值多项式 ,用于离散数据的拟合 C/C++ code#include #include #include float lagrange(float *x,float *y,float xx,int n) /*拉格朗日插值算法 */ { int i,j; float *a,yy=0.0; /*a作为临时变量,记录拉格朗日插值多项式 */ a=(float *)malloc(n*sizeof(float)); for(i=0;i<> { a[i]=y[i]; for(j=0;j<> if(j!=i) a[i]*=(xx-x[j])/(x[i]-x[j]); yy+=a[i]; } free(a); return yy; } main() { int i,n; float x[20],y[20],xx,yy; printf( scanf( if(n>=20) {printf( { printf( scanf( } printf( for(i=0;i<> { printf( printf( printf( scanf( yy=lagrange(x,y,xx,n); printf( getch(); } 2.牛顿插值多项式,用于离散数据的拟合 C/C++ code#include #include #include void difference(float *x,float *y,int n) { float *f; int k,i; f=(float *)malloc(n*sizeof(float)); for(k=1;k<> { f[0]=y[k]; for(i=0;i<> f[i+1]=(f[i]-y[i])/(x[k]-x[i]); y[k]=f[k]; } return; } main() { int i,n; float x[20],y[20],xx,yy; printf( scanf( if(n>=20) {printf( { printf( scanf( } printf( for(i=0;i<> { printf( printf( difference(x,(float *)y,n); printf( scanf( yy=y[20]; for(i=n-1;i>=0;i--) yy=yy*(xx-x[i])+y[i]; printf( getch(); } 3.高斯列主元消去法 , 求解其次线性方程组 C/C++ code#include #include #define N 20 int main() { int n,i,j,k; int mi,tmp,mx; float a[N][N],b[N],x[N]; printf( scanf( if(n>N) { printf( getch(); return 1; } if(n<> { printf( getch(); return 1; } printf( { for(j=0;j<> scanf( printf( scanf( for(i=0;i<> { for(j=i+1,mi=i,mx=fabs(a[i][j]);j { mi=j; mx=fabs(a[j][i]); } if(i<> { tmp=b[i];b[i]=b[mi];b[mi]=tmp; for(j=i;j<> { tmp=a[i][j]; a[i][j]=a[mi][j]; a[mi][j]=tmp; } } for(j=i+1;j<> { tmp=-a[j][i]/a[i][i]; b[j]+=b[i]*tmp; for(k=i;k<> a[j][k]+=a[i][k]*tmp; } } x[n-1]=b[n-1]/a[n-1][n-1]; for(i=n-2;i>=0;i--) { x[i]=b[i]; for(j=i+1;j<> x[i]-=a[i][j]*x[j]; x[i]/=a[i][i]; } for(i=0;i<> printf( return 0; } #include #include #define NUMBER 20 #define Esc 0x1b #define Enter 0x0d float A[NUMBER][NUMBER+1] ,ark; int flag,n; exchange(int r,int k); float max(int k); message(); main() { float x[NUMBER]; int r,k,i,j; char celect; clrscr(); printf( printf( printf( celect=getch(); if(celect==Esc) exit(0); printf( scanf( printf( for(i=1;i<> { printf( for(j=1;j<=n+1;j++)> for(k=1;k<> { ark=max(k); if(ark==0) { printf( } else if(flag!=k) exchange(flag,k); for(i=k+1;i<> for(j=k+1;j<> A[i][j]=A[i][j]-A[k][j]*A[i][k]/A[k][k]; } x[n]=A[n][n+1]/A[n][n]; for( k=n-1;k>=1;k--) { float me=0; for(j=k+1;j<> { me=me+A[k][j]*x[j]; } x[k]=(A[k][n+1]-me)/A[k][k]; } for(i=1;i<> { printf( } message(); } exchange(int r,int k) { int i; for(i=1;i<> A[0][i]=A[r][i]; for(i=1;i<> A[r][i]=A[k][i]; for(i=1;i<> A[k][i]=A[0][i]; } float max(int k) { int i; float temp=0; for(i=k;i<> if(fabs(A[i][k])>temp) { temp=fabs(A[i][k]); flag=i; } return temp; } message() { printf( { case Enter: main(); case Esc: exit(0); default:{printf( } 4.龙贝格求积公式,求解定积分 C/C++ code#include #include #define f(x) (sin(x)/x) #define N 20 #define MAX 20 #define a 2 #define b 4 #define e 0.00001 float LBG(float p,float q,int n) { int i; float sum=0,h=(q-p)/n; for (i=1;i<> sum+=f(p+i*h); sum+=(f(p)+f(q))/2; return(h*sum); } void main() { int i; int n=N,m=0; float T[MAX+1][2]; T[0][1]=LBG(a,b,n); n*=2; for(m=1;m<> { for(i=0;i<> T[i][0]=T[i][1]; T[0][1]=LBG(a,b,n); n*=2; for(i=1;i<> T[i][1]=T[i-1][1]+(T[i-1][1]-T[i-1][0])/(pow(2,2*m)-1); if((T[m-1][1] return ; } } } C/C++ code5.牛顿迭代公式,求方程的近似解 C/C++ code#include #include #include #define N 100 #define PS 1e-5 #define TA 1e-5 float Newton(float (*f)(float),float(*f1)(float),float x0 ) { float x1,d=0; int k=0; do { x1= x0-f(x0)/f1(x0); if((k++>N)||(fabs(f1(x1))<> { printf( getch(); exit(); } d=(fabs(x1)<1?x1-x0:(x1-x0)> x0=x1; printf( } while((fabs(d))>PS&&fabs(f(x1))>TA) ; return x1; } float f(float x) { return x*x*x+x*x-3*x-3; } float f1(float x) { return 3.0*x*x+2*x-3; } void main() { float f(float); float f1(float); float x0,y0; printf( scanf( printf( y0=Newton(f,f1,x0); printf( getch(); } 6. 牛顿 -科特斯求积公式,求定积分 C/C++ code#include #include int NC(a,h,n,r,f) float (*a)[]; float h; int n,f; float *r; { int nn,i; float ds; if(n>1000||n<> { if (f) printf( return(-1); } if(n==2) { *r=0.5*((*a)[0]+(*a)[1])*(h); return(0); } if (n-4==0) { *r=0; *r=*r+0.375*(h)*((*a)[n-4]+3*(*a)[n-3]+3*(*a)[n-2]+(*a)[n-1]); return(0); } if(n/2-(n-1)/2<> nn=n; else nn=n-3; ds=(*a)[0]-(*a)[nn-1]; for(i=2;i<> ds=ds+4*(*a)[i-1]+2*(*a)[i]; *r=ds*(h)/3; if(n>nn) *r=*r+0.375*(h)*((*a)[n-4]+3*(*a)[n-3]+3*(*a)[n-2]+(*a)[n-1]); return(0); } main() { float h,r; int n,ntf,f; int i; float a[16]; printf( for(i=0;i<> scanf( h=0.2; f=0; ntf=NC(a,h,n,&r,f); if(ntf==0) printf( else printf( getch(); } 7.雅克比迭代,求解方程近似解 C/C++ code#include #include #define N 20 #define MAX 100 #define e 0.00001 int main() { int n; int i,j,k; float t; float a[N][N],b[N][N],c[N],g[N],x[N],h[N]; printf( if(n>N) { printf( {printf( for(i=0;i<> for(j=0;j<> scanf( printf( for(i=0;i<> scanf( for(i=0;i<> for(j=0;j<> { b[i][j]=-a[i][j]/a[i][i]; g[i]=c[i]/a[i][i]; } for(i=0;i<> { for(j=0;j<> h[j]=g[j]; { for(k=0;k<> { if(j==k) continue; h[j]+=b[j][k]*x[k]; } } t=0; for(j=0;j<> if(t for(j=0;j<> x[j]=h[j]; if(t<> { printf( for(i=0;i<> printf( getch(); return 0; } printf( return 1; } getch(); } 8.秦九昭算法 C/C++ code#include float qin(float a[],int n,float x) { float r=0; int i; for(i=n;i>=0;i--) r=r*x+a[i]; return r; } main() { float a[50],x,r=0; int n,i; do { printf( scanf( } while(n<> printf( for(i=0;i<> scanf( printf( scanf( r=qin(a,n,x); printf( getch(); } 9.幂法 C/C++ code#include #include #define N 100 #define e 0.00001 #define n 3 float x[n]={0,0,1}; float a[n][n]={{2,3,2},{10,3,4},{3,6,1}}; float y[n]; main() { int i,j,k; float xm,oxm; oxm=0; for(k=0;k<> { for(j=0;j<> { y[j]=0; for(i=0;i<> y[j]+=a[j][i]*x[i]; } xm=0; for(j=0;j<> if(fabs(y[j])>xm) xm=fabs(y[j]); for(j=0;j<> y[j]/=xm; for(j=0;j<> x[j]=y[j]; if(fabs(xm-oxm)<> { printf( printf( for(k=0;k } oxm=xm; } getch(); } 10.高斯塞德尔 C/C++ code#include #include #define N 20 #define M 99 float a[N][N]; float b[N]; int main() { int i,j,k,n; float sum,no,d,s,x[N]; printf( scanf( if(n>N) { printf( return 1; } if(n<> { printf( for(i=0;i<> for(j=0;j<> scanf( printf( for(i=0;i for(i=0;i k=0; printf( for(i=0;i do { k++; if(k>M){printf( break; } no=0.0; for(i=0;i<> { s=x[i]; sum=0.0; for(j=0;j<> if (j!=i) sum=sum+a[i][j]*x[j]; x[i]=(b[i]-sum)/a[i][i]; d=fabs(x[i]-s); if (no } printf( for(i=0;i } while (no>=0.1e-6); if(no<> { printf( printf( for (i=0;i<> printf( } getch(); } 《MATLAB程序设计实践》 1、编程实现以下科学计算算法,并举一例应用之。 “中点公式法和五点公式法求数值微分” 解:中点公式法和五点公式法求数值微分如下: 例5-4:中点公式法求导数应用实例。采用中点公式法求函数 xf=在x=4处的导数。 解:在MATLAB命令窗口中输入: >>df=MidPoint('sqrt(x)',4) 输出结果为: df= 0.2500 x采用中点公式法求函数f=在x=4处的导数为0.25，而导数的精确值也是0.25 . 详见以下, 中点公式法流程图, 开始 If nargin=2,h=0.1 if nargin=3,h=0判断输入参 数个数 输出:h用中点公式 结束 不能为求数值微分 0 源代码: function df=MidPoint(func,x0,h) if nargin == 2 h = 0.1; else if (nargin == 3 && h == 0.0) disp('h???ü?a0??'); return; end end y1 = subs(sym(func), findsym(sym(func)),x0+h); y2 = subs(sym(func), findsym(sym(func)),x0-h); df = (y1-y2)/(2*h); 运行结果如下, 例5-5:五点公式法求导数应用实例。采用五点公式法求函数 f=sin(x)在x=2处的导数。 解:在MATLAB命令窗口中输入: >>df1=FivePoint('sin(x)',2,1); >>df2=FivePoint('sin(x)',2,2); >>df3=FivePoint('sin(x)',2,3); >>df4=FivePoint('sin(x)',2,4); >>df5=FivePoint('sin(x)',2,5); 用五种方法得到的结果为: df1= -0.4161 df2= -0.4161 df3= -0.4161 df4= -0.4161 df5= -0.4161 而函数在f=sin(x)在x=2的导数为cos(2)=-0.4161,从上面的结果来看，五点公式的 精度是很高的。 详见以下, 五点公式法流程图, 开始 判断输入参 If nargin=3,h=0.1 if nargin =4,h=0 数个数 输出: h 不能为用公用公用公用公用公 0 式一式二式三式四式五 求导 求导 求导 求导 求导 结束 源代码: function df=FivePoint(func,x0,type,h) %??μ???ê???,?óè?o?êyfunc?úx0??μ?êy ×?íòê?è?òa %o?êy???ofunc %?óμ?μ??ox0 %??ê?μ?D?ê??otype??è?1,2,3,4,5,?? %à?é??????oh %μ?êy?μ?odf if nargin ==3 h=0.1; else if (nargin ==4&&h==0.0) disp('h???ü?a0'); return; end end y0 = subs(sym(func),findsym(sym(func)),x0); y1 = subs(sym(func),findsym(sym(func)),x0+h); y2 = subs(sym(func),findsym(sym(func)),x0+2*h); y3 = subs(sym(func),findsym(sym(func)),x0+3*h); y4 = subs(sym(func),findsym(sym(func)),x0+4*h); y_1 = subs(sym(func),findsym(sym(func)),x0-h); y_2 = subs(sym(func),findsym(sym(func)),x0-2*h); y_3 = subs(sym(func),findsym(sym(func)),x0-3*h); y_4 = subs(sym(func),findsym(sym(func)),x0-4*h); switch type case 1, df=(-25*y0+48*y1-36*y2+16*y3-3*y4)/(12*h);%ó?μúò?????ê??óμ?êy case 2, df=(-3*y_1-10*y0+18*y1-6*y2+y3)/(12*h);%ó?μú?t????ê??óμ?êy case 3, df=(y_2-8*y_1+8*y1-y2)/(12*h);%ó?μúèy????ê??óμ?êy case 4, df=(3*y1+10*y0-18*y_1+6*y_2-y_3)/(12*h);%ó?μú??????ê??óμ?êy case 5, df=(25*y0-48*y_1+36*y_2-16*y_3+3*y_4)/(12*h);%ó?μú??????ê??óμ?êy end 运行结果如下, 2、编程解决以下科学计算和工程实际问题。 ?已知阿波罗(Apollo)卫星的运动轨迹(x,y)满足下列微分方程 **..，，,x，,,(x,,)， x,2y，x,,33 rr12 *...,y,y,,2，,, yxy33 rr12 22*22*1,(x，,)，y,(x，,)，y其中= ，=1- ， 试在初值,,,rr1282.45 .. x(0)=1.2, , 下进行数值求解，并绘制出阿波罗卫星位置(x,y)x(0),0y(0),,1.04935751, 的轨迹。 ?解:根据题目选用MATLAB代码如下: function dy=weifen(t,y) % 编程解决阿波罗(Apollo)卫星的运动轨迹 求解器属于变步长的一种，采用Runge-Kutta算法 万事如意 u=1/82.45; b=1-u; dy=zeros(4,1); r=zeros(2,1); r(1)=sqrt((y(1)+u)^2+(y(3))^2); r(2)=sqrt((y(1)+b)^2+(y(3))^2); dy(1)=y(2); dy(2)=2*dy(3)+y(1)-b*(y(1)+u)/(r(1)^3)-u*(y(1)-b)/(r(2)^3); dy(3)=y(4); dy(4)=-2*dy(1)+y(3)-b*y(3)/(r(1)^3)-u*y(3)/(r(2)^3); 解:在MATLAB命令窗口中输入 >>ode45('weifen',[0 2.00],[1.2 0 0 -1.04935751]) >>[T,Y]=ode45('weifen',[0 1.26],[1.2 0 0 -1.04935751]) 运行结果: 阿波罗卫星位置(x,y)的轨迹图如下: ?实验图所示是一个跷跷板，两板价较为，左边板长为1.5m，上面的小孩重150N,右边板长为2m,小孩重为400N.求当跷跷板平衡时，左边木板与水平方向的夹角ɑ的大小。要求先求解析解，然后给出两种解决方案。 ?解:根据力矩平衡求解析解 由图示可有下列关系式: 1cos(,,,)cos,5001.5=2400 ，，3 解该式得: ,,350cos,4003sin 83, tan,7 83,,arctan7 ,,0.4678rad即: 两种方法的求解: 方案一:采用两步迭代法求解方程: 1cos(,,,) 5001.5=2400 cos,，，3 两步迭代法的MATLAB的代码如下: 源代码, function root=TwoStep(f,a,b,type,eps) if(nargin==4) eps=1.0e-4; %eps ?íê?μü?ú???è end f1=subs(sym(f),findsym(sym(f)),a); f2=subs(sym(f),findsym(sym(f)),b); if(f1==0) root=a; end if(f2==0) root=b; end if(f1*f2>0) disp('两端点函数值乘积大于0!'); return; else tol=1; fun=diff(sym(f)); fa=subs(sym(f),findsym(sym(f)),a); fb=subs(sym(f),findsym(sym(f)),b); dfa=subs(sym(fun),findsym(sym(fun)),a); dfb=subs(sym(fun),findsym(sym(fun)),b); if(dfa>dfb) root=a; else root=b; end while(tol>eps) if(type==1) r1=root; fx1=subs(sym(f),findsym(sym(f)),r1); dfx=subs(sym(fun),findsym(sym(fun)),r1); r2=r1-fx1/dfx; fx2=subs(sym(f),findsym(sym(f)),r2); root=r2-fx2/dfx; tol=abs(root-r1); else r1=root; fx1=subs(sym(f),findsym(sym(f)),r1); dfx=subs(sym(fun),findsym(sym(fun)),r1); r2=r1-fx1/dfx; fx2=subs(sym(f),findsym(sym(f)),r2); root=r2-fx2*(r2-r1)/(2*fx2-fx1); tol=abs(root-r1); end end end 解:在MATLAB命令窗口中输入 r=TwoStep('500*1.5*cos(x)-2*400*cos(1/3*pi-x)',0,1/3*pi,1) 运行结果如下: 两个小孩所产生力矩随变化的曲线如下图片: , 运用Data cursor 工具,得到交点处对应的X值为:0.4678 ,,0.4678也即, >> x=0:0.0001:1/3*pi; >> ML=500*1.5*cos(x); >> MR=400*2*cos(1/3*pi-x); ,'-',x,MR ,'-') >> plot(x,ML 中点坐标公?式在直线方?程中的应用 ? 陕西汉中市?405学校? 侯有岐 72331?2 直线方程问?题是学习解?析几何的基?础,涉及到的知?识比较多,若能灵活应?用,将收到事半?功倍的效果?.本文例谈中?点坐标公式?在直线方程?中的应用,仅供参考. 一、 过一点直线?夹在两条已?知直线间的?线段中点问?题 l例1:直线 被两条直线? 和 截得的线l: 4x，y，3,0l: 3x,5y,5,012 l段?中点为P，求直线 的方程. (,1,2) l分析: 本题同学们?很可能设出? 的斜截式,由方程组解?得 的交点 l 与 l、l12 ,5,kk，85k，158k，6k,,3A(,)、 B(，), 再由中点公?式解得 ,代入斜截式?k，4k，43,5k3,5k 即可,但计算过于?繁琐.我们可利用?中点坐标,设出A、B两点坐标?来求. 解: 设 , l 与 l 交于A(a，b)，l 与 l 交于B12 因为P是 AB 中点,所以B点为? , (,2,a,4,b) 4a，b，3,0a,,2,, 于是, 解得 即 A(,4， 2),,3(,2,a),5(4,b),5,0b,5,, l 由两点式知?所求直线 的方程为 . 3x，y，1,0 点评: 由本例可知?,求出的交点?,虽然思路简?单,但计算过于?繁琐.面l 与 l、l12 对此类问?题,先设出再利?用中点坐标?公式表示出?l 与 l 交于A(a，b)，1 ,由 A、B两点分别?在上，l 与 l 交于B(2m,a,2n,b)(设中点P(m，n))l、l212 (a，b)通过解方程?组求出点A?，后由直线的?两点式或点?斜式即可求?得所求直线?方程. 二、 利用中点坐?标公式求解?决有关对称?问题 例2:(1)已知点A(5，8),,求A点关于?B点的对称?点C的坐标?. B(4,1) ll (2) 已知点A(,直线 方程为 , 求点A关于?的对,4,4)3x，y,2,0 ,A称点的?坐标. l,1 (3) 求 3x,y,4,0 关于点P (2,)对称的直线? 的方程. 分析: (1)可直接应用?中点坐标公?式解决; (2)可根据点关?于线对称的?特点,利用垂直平?分解方程组?解决; (3)可利用中点?坐标公式,由中心对称?的定义解决?. 5，x,4,,x,3,,2解: (1) 设C(, 由中点坐标?公式有 得 x,y),,8，yy,,6,,1,,2, 所以点C为?所求. (3,,6) ,,,,,,lAA,l (2) 设点,因为点A与?关于 对称,所以 , 且中点AAAA(x,y) ,y,4,,(,3),,1,,x,2,,,x，4l在 上,所以 解得 ,,,,,x,4y，4y,6,,3,，,2,0,22, , 所以点为所?求. A(2,6) l (3) 设 上任意一点?为,关于P(2,1) 对称点()在直线(x,y)4,x,,2,y 上, 所以 , 3x,y,4,03(4,x),(,2,y),4,0 l 所以 , 则所求直线? 为 . 3x,y,10,03x,y,10,0 点评: 涉及“点关于点对?称”，“点关于线对?称”，“线关于点对?称”，“线关于线对?称”等问题时，若能恰当应?用中点坐标?公式，就能使问题?的解决更加?简洁和富有?创新。 ,l 一般情况下?，点E关于直?线 : 的对称点，E(x,y)(a,b)Ax，By，C,000可通过垂直?平分的特点?，利用中点坐?标公式，解方程组 ,ybA,0,(,),,1,,xaB,0 解决. ,，，xayb,00,，,，,0ABC,22, 借助中点坐标公式，理解对称性的性质 贵州省天柱民族中学 杨武灵 （电话：13985293350 邮编：556600 E-mail：tzmz150@sohu.com） [摘要] 借助中点坐标公式，帮助理解函数的对称性性质，使对称性 问题更简捷地得到解决，有利于提高学生的学习兴趣。 [关键词] 中点坐标公式 函数 对称性 我们知道，函数是中学数学教学中的重要内容，也是进一步学习与应用的重要基础，更是高考的重点与亮点，而对称性却是函数的一个基本性质，也是教学中的一个难点。在教材中与对称性有关的问题从函数的奇偶性、周期性中得以体现，并且广泛存在于数学问题之中，充分体现了数学之美，利用好对称性，往往能使问题更简捷地得到解决，有利于提高学生的学习兴趣。尽管我们总结了许多条对称性的性质，也很难使学生突破这一难点。本文借助中点坐标公式，很容易理解与对称有关的图形性质。 一、对称性的理解 在理解对称性的性质时，先应理解点的坐标意义，对于点（x ，b ）中，b 是常数，x 是变量，故点（x ，b ）的图像是直线y= b；同理点（a ，y ）的图像是直线x= a ；而（a ，b ）为一定点。这样，借助中点坐标公式，求出A 与B 的中点坐标，很容易得到以下性质。 1、性质1：点A （x ，y ）与B （x ，2b －y ）关于直线y= b对称. 推论1：点（x ，y ）关于x 轴（直线y=0）对称的点为（x ，－y ）. 推论2：函数 y = f (x)与y =－f (x)的图像关于x 轴（直线y=0） 对称. 2、性质2：点A （x ，y ）与B （2a －x ，y ）关于直线x= a对称. 推论1：点（x ，y ）关于y 轴（直线x=0）对称的点为（－x ， y ）. 推论2：函数 y = f (x)的图像关于y 轴（直线x=0）对称的充要条件是f (x) = f (－x). 推论3：函数 y = f (x)的图像关于直线x = a对称的充要条件是 f (a+x) = f (a－x) 或f (x) = f (2a－x). 3、性质3：点A(x，y) 与B(2a－x ，2b －y) 关于点（a ，b ）对称。 推论1：点（x ，y ）关于原点O 对称的点为（－x ，－y ）. 推论2：函数 y = f (x)的图像关于原点O 对称的充要条件是f ( x) +f (－x) =0. 推论3：函数 y = f (x)的图像关于（a ，b ）对称的充要条件是f (2a－x) =2b－f (x). 4、性质4：若点A （x ，y ）关于直线y =ax+b对称的点为B （x 0，y 0），则由A 、B 的中点在直线y =ax+b上，且AB 与直线y =ax+b垂直得：y+ y0=a（ x+ x0）+2b，x －x 0=－a （y －y 0），从而解得： x 0=-2a （ax -y +b ）（2ax -y +b ）+x +y . ， y =0a 2+1a 2+1 推论1：点（x ，y ）关于直线y = x对称的点为（y ，x ）. 推论2：点（x ，y ）关于直线y =－ x 对称的点为（－y ，－x ）. 推论3：函数 y = f (x)与其反函数y = f－1 (x)（即x = f (y)）的图像关于直线y = x对称 5、对于曲线f (x，y)=0，若A （x ，y ）关于某点（或直线）的 对称点为B （x 0，y 0），则有： （1）曲线 f (x，y)=0 关于该点（或直线）对称的曲线方程是 f (x0，y 0) = 0. （2）曲线 f (x，y)=0 关于该点（或直线）对称的充要条件是 f (x，y)= f (x0，y 0) . 二、 函数对称性应用举例 例1、定义在R 上的非常数函数满足：f (5+x)为偶函数，且f (10－x) = f (10+x)，则f (x)一定是（ ） (A)是偶函数，也是周期函数 (B)是偶函数，但不是周期函数 (C)是奇函数，也是周期函数 (D)是奇函数，但不是周期函数 解：由题意知，f (5+x) 向右平移5个单位后得f (x)，所以x = 5是f (x)的对称轴，即f (x)= f (10－x) ，又因为f (10－x)=f (10+x)，所以x = 10也是f (x)的对称轴，且f(x) = f (10－x)= f (10+x) ，因此f (x)是以10为其一个周期的周期函数， 从而x =0即y 轴也是f (x)的对称轴。 故选(A) 例2、设f(x)是定义在R 上的奇函数，且f(2+x)=－f(x)，当－1≤x ≤0时，f (x) =2x，则f (8.5) = ____. 解：依题意知，f(2+x)=－f(x) = f(－x) ，所以f(4+x)=－f(2+x) = f( x)， 故y = f(x)是以4为周期的周期函数，可得：f (8.5 ) = f (8+0.5 ) = f (0.5 ) = －f (－0.5 ) =1. 例3. 以x+2y+1=0为对称轴，直线x －y －2=0的轴对称图形的方程是（ ） (A)7x+y－8=0 (B) 7x－y+8=0 (C) 7x－y －8=0 (D) 7x+y+8=0 解：由性质4易知，（x ，y ）关于x+2y+1=0的对称点为 P （x -y -，-x -y -），将P 点坐标代入x －y －2=0后化 简得：7x －y －8=0，故选(C). 例4. （2004年安徽）已知抛物线C 1：y=ax2+bx与C 2：y= 2px （p ＞0）关于x+y=1对称，（1）求a 、b 、p 的值；（2）求焦点C 1与C 2间的距离. 解：（1）法1：依题意得：C 2的顶点为（0，0），可求得C 2关于x+y=1的对称点（1，1）为C 1的顶点，故点（1，1）同时在抛物线 ?a +b =1?b 1C 1、C 2上，所以有：?=1 从而解得：a=－1，b=2，p=. ?-2?2a ??1=2p 354525453545 法2：由于（x ，y ）关于x+y=1的对称点为（1－y ，1－x ），所以（1－x ）2=2p（1－y ）与y=ax2+bx是相同的曲线，比较两曲线的系数得a=－1，b=2，p=. （2）由于C 2的焦点坐标为（，0），它关于x+y=1的对称点为 22（1 ， ），所以两焦点的距离为1-）+（-0）=12143 4143432. 4 例5、证明：若函数y = f (x) 图像同时关于点A (a1 ,b)和B (a2 ,b)成中心对称（a 1≠a 2），则y = f (x)是周期函数，其中一个周期为2(a1－a 2). 证明： 由性质3得：f (2a1－x) =2b－f (x)，f (2a2－x) =2b－f (x). 所以f (2a1－x) =f (2a2－x) ， 从而f (x +2a1－2a 2) =f (2a1－(2a2－x)) =f (2a2－(2a2－x)) = f (x)，即y = f (x)是周期函数，其中一个周期为2( a1－a 2). 例6、已知：函数y = f (x)关于 点(1 ,4) 成中心对称，又关于直线x = 3成轴对称，当1≤x ≤2时，f (x) = x，求f (2009)的值. 解：由题意知，f (2－x) =8－f (x)，f (6－x) = f (x)，则有 f (x) =f (2－(2－x )) = 8－f (2－x) =8－f (6－(4+x)）=8－f(4+x) =8－f( 2－(－2－x )) = f (－2－x) =f (6－(8+x)）=f (8+x). 即周期T=8. 故f (2009) =f (8×251+1) =f (1) = 1. 转载请注明出处范文大全网 » 中点公式法和五点公式法求数值微分中点公式法和五点公式法求数值微分[终稿]
中点坐标公式在直线方程中的应用
借助中点坐标公式
