本人就读于某211大学数学系，本着造福学弟学妹的原则上传了自己的代码，仅供参考
原本是专业必修课用来交作业的c语言版本，但是小学期老师要MATLAB版的，看着身边同学一个个抓耳挠腮，所以更新了MATLAB版本的代码，适用于那些不想申请A的学弟学妹。（申请A老师会细看你的代码，会发现是抄的）
同校的自然一看就懂我在说什么
牛顿插值法（课本32页例4）
c语言版本
#include <stdio.h>
#include <stdlib.h>
int main()
{
float table(int n,float a1[10],float a2[10],float a3[10][10]);
float newton(int n,float a4[10][10],float a5[10]);
float arrX[10],arrY[10],arrL[10][10];
int num,i;
printf("请输入插值点个数（个数应小于10）");
scanf("%d",&num);
printf("请输入各个插值节点的值:\n");
for(i=0; i<num; i++)
{
printf("请输入X%d值:",i+1);
scanf("%f",&arrX[i]);
printf("请输入Y%d值:",i+1);
scanf("%f",&arrY[i]);
}
table(num,arrX,arrY,arrL);
newton(num,arrL,arrX);
return 0;
}
float table(int n,float a1[10],float a2[10],float a3[10][10])
{
int i,j;
for(i=0; i<n; i++)
{
a3[i][0]=a2[i];//第一列初始化为y值
}
for(j=0; j<n; j++)
for(i=n-1; i>j; i--)//从第一列往前循环，到i=j为止，即下三角全部计算赋值
{
a3[i][j+1]=(a3[i][j]-a3[i-1][j])/(a1[i]-a1[i-1*(j+1)]);//差商表计算，最后一项arrX[i-1*(j+1)]角标是计算步长
}
for(j=1; j<n; j++)//上三角未计算赋值，系统会随机分配值给上三角的每一位，故赋值0会使差商表输出更整齐
for(i=0; i<j; i++)
{
a3[i][j]=0.0;
}
printf("差商表为:\n");
printf("\n");
for(i=0; i<n; i++)
{
for(j=0; j<n; j++)
{
if(j%n==0)//num个数一列输出
printf("\n");
printf("%f\t",a3[i][j]);
}
}
printf("\n");
printf("\n");
return 0;
}
float newton(int n,float a4[10][10],float a5[10])
{
int i;
float x,y,t1=1.0;
while(1)
{
printf("请输入要插入节点的x值:");
scanf("%f",&x);
y=a4[0][0];
for(i=1; i<n; i++)
{
t1=t1*(x-a5[i-1]);
y=y+a4[i][i]*t1;//差商表对角线值依次乘以(x-x0)(x-x1)。。。。。。
}
printf("插值结果为:%f",y);
printf("\n");
}
return 0;
}
MATLAB版本(data.txt为同级目录新建文本文件，把题目已知数据输入，注意空格和换行)
close all;clc;clear;
fileID=fopen('data.txt');
A=textscan(fileID,'%s %s');
fclose(fileID);
x=A{1,1};
y=A{1,2};
format long
x=str2double(x);
y=str2double(y);
a=[x,y];
[m,n]=size(a);
a=([a zeros(m,m-n)])+zeros(m,m);
for i=2:1:m
for j=3:1:i+1
a(i,j)=(a(i,j-1)-a(i-1,j-1))/(a(i,1)-a(i-j+2,1));
end
end
a
x=input('请输入x：','s');
x=str2double(x);
n=input('请输入插值多项式的次数n：','s');
d=a(1,2);
for i=2:1:m
b=1;
for j=3:1:i+1
b=(x-a(i-j+2,1))*b;
c=a(i,j)*b;
end
d=d+c;
end
y=sprintf('近似的结果f(x)= %.5f',d);
disp(y)
数值积分章节的复合梯形和复合辛普森公式（课本108页例3）
c语言版本
#include "stdio.h"
#include "math.h"
void main()
{
int k;
double a,b,n1,n2,h,x;
double f,f1,f2,T,F=0.0;
printf("请输入积分区间a b:");
scanf("%lf%lf",&a,&b);
printf("请输入使用复合梯形公式时将区间划分的份数n1:");
scanf("%lf",&n1);
printf("请输入使用复合辛普森公式时将区间划分的份数n2:");
scanf("%lf",&n2);
h=(b-a)/n1;
for(k=1;k<n1;k++)
{
x=a+k*h;
f=sin(x)/x;
F+=f;
}
if(a!=0.0)
f1=sin(a)/a;
else(f1=1.0);
f2=sin(b)/b;
T=0.5*h*(f1+2*F+f2);
printf("\n使用复合梯形公式求得T=%0.7f",T);
double F1=0.0,F2=0.0,S=0.0;
h=(b-a)/n2;
for(k=0;k<n2;k++)
{
x=a+k*h+0.5*h;
f=sin(x)/x;
F1+=f;
}
for(k=1;k<n2;k++)
{
f=0;
x=a+k*h;
f=sin(x)/x;
F2+=f;
}
if(a!=0)
f1=sin(a)/a;
else
f1=1.0;
f2=sin(b)/b;
S=h/6.0*(f1+4*F1+2*F2+f2);
printf("\n使用复合辛普森公式求得S=%0.7f",S);
}
MATLAB版本
close all;clc;clear;
a=input('请输入积分下限a：','s');
a=str2double(a);
b=input('请输入积分上限b：','s');
b=str2double(b);
n1=input('请输入使用复合梯形公式时将区间划分的份数n1：','s');
n1=str2double(n1);
n2=input('请输入使用复合辛普森公式时将区间划分的份数n2：','s');
n2=str2double(n2);
F=0;
h=(b-a)/n1;
for k=1:1:n1-1
x=a+k*h;
f=sin(x)/x;
F=F+f;
end
if a~=0.0
f1=sin(a)/a;
else
f1=1.0;
end
f2=sin(b)/b;
T=0.5*h*(f1+2*F+f2);
y=sprintf('使用复合梯形公式求得T= %.7f',T);
disp(y)
F1=0.0;
F2=0.0;
S=0.0;
h=(b-a)/n2;
for k=0:1:n2-1
x=a+k*h+0.5*h;
f=sin(x)/x;
F1=F1+f;
end
for k=1:1:n2-1
f=0;
x=a+k*h;
f=sin(x)/x;
F2=F2+f;
end
if a~=0
f1=sin(a)/a;
else
f1=1.0;
end
f2=sin(b)/b;
S=h/6.0*(f1+4*F1+2*F2+f2);
y=sprintf('使用复合辛普森公式求得S= %.7f',S);
disp(y)
线性方程组的直接解法之追赶法（书177页习题9）
c语言版本
#include "stdio.h"
#include "math.h"
#define n 5
void main()
{
float a[n]={0,-1,-1,-1,-1},b[n]={2,2,2,2,2},c[n]={-1,-1,-1,-1,0},
s[n],t[n],temp=0,x[n]={1,0,0,0,0};
int k;
for(k=0;k<n;k++)
{
s[k]=b[k]-a[k]*temp;
t[k]=c[k]/s[k];temp=t[k];
}
temp=0;
for(k=0;k<n;k++)
{
x[k]=(x[k]-a[k]*temp)/s[k];
temp=x[k];
}
for(k=n-2;k>=0;k--)
x[k]=x[k]-t[k]*x[k+1];
printf("\nx= \n");
for(k=0;k<n;k++)
printf(" %f \n",x[k]);
}
MATLAB版本（data1.xlsx、data2.xlsx为A和b，下同，其实这是正确导入数据的方式，第一个程序里用txt文件存数据不规范）
close all;clc;clear;
A = xlsread('data1.xlsx');
x = xlsread('data2.xlsx');
[m,n]=size(A);
a=[0;diag(A,-1)];
b=diag(A,0);
c=[diag(A,1);0];
temp=0;
for k=1:1:m
s(k,1)=b(k,1)-a(k,1)*temp;
t(k,1)=c(k,1)/s(k,1);
temp=t(k,1);
end
temp=0;
for k=1:1:m
x(k,1)=(x(k,1)-a(k,1)*temp)/s(k,1);
temp=x(k,1);
end
for k=n-1:-1:1
x(k,1)=x(k,1)-t(k,1)*x(k+1,1);
end
sprintf(' %.6f',x)
线性方程组的迭代法之高斯-赛德尔迭代法（课本189页例6）
c语言版本
#include "stdio.h"
#include "math.h"
#define n 3
void main()
{
float A[n][n]={{8,-3,2},{4,11,-1},{6,3,12}},
b[n]={20,33,36},x[n]={0.0,0.0,0.0},epsilon,normal,temp,t;
int N,i,j,k;
printf("\n 精度 = "); scanf("%f",&epsilon);
printf("\n 最大迭代精度N = "); scanf("%d",&N);
printf("\n %d : ",0);
for(i=0;i<n;i++)
printf("%f",x[i]);
for(k=0;k<N;k++)
{
normal=0;
for(i=0;i<n;i++)
{
t=x[i];
x[i]=b[i];
for(j=0;j<n;j++)
if(j!=i)
x[i]=x[i]-A[i][j]*x[j];
x[i]=x[i]/A[i][i];
temp=fabs(x[i]-t);
if(temp>normal)
normal=temp;
}
printf("\n %d :",k+1);
for(i=0;i<n;i++)
printf(" %f",x[i]);
if(normal<epsilon)
{
printf("\n");
return;
}
}
printf("\n \n迭代%d次后仍未求得满足精度的解\n",N);
}
MATLAB版本
close all;clc;clear;
A = xlsread('data1.xlsx');
b = xlsread('data2.xlsx');
[m,n]=size(A);
x=[0;0;0];
e=input('请输入精度e：','s');
e=str2double(e);
N=input('请输入最大迭代精度N：','s');
N=str2double(N);
for i=1:1:m
sprintf('%f',x(i,1))
end
for k=1:1:N
normal=0;
for i=1:1:m
t=x(i,1);
x(i,1)=b(i,1);
for j=1:1:m
if j~=i
x(i,1)=x(i,1)-A(i,j)*x(j,1);
end
end
x(i,1)=x(i,1)/A(i,i);
temp=abs(x(i,1)-t);
if(temp>normal)
normal=temp;
end
end
for i=1:1:m
sprintf('%f',x(i,1))
end
if normal<e
return;
end
end
常微分初值问题三种数值解法(课本284页例2三种解法)
c语言版本
#include "stdio.h"
#include "math.h"
//
#define f(x,y) ((y)-2*(x)/(y))
#define y(x) sqrt(1+2*(x))
void main()
{
float a=0,b=1.0f,h=0.1f,y0=1.0f,x,ye,yp,ym,k1,k2,k3,k4,yr,yx;
printf("\n 欧拉方法");
printf(" 改进欧拉方法 四阶龙格库塔方法");
printf("\nx ye[k] ");
printf("ym[k] yr[k] \n");
printf("%3.1f %9.6f ",a,0);
printf("%9.6f %9.6f \n",y0,y0);
x=a;
ye=y0;
ym=y0;
yr=y0;
while(x<b)
{
ye=ye+h*f(x,ye);
yp=ym+h*f(x,ym);
ym=ym+h/2*(f(x,ym)+f(x+h,yp));
k1=f(x,yr);
k2=f(x+h/2,yr+h/2*k1);
k3=f(x+h/2,yr+h/2*k2);
k4=f(x+h,yr+h*k3);
yr=yr+h/6*(k1+2*k2+2*k3+k4);
x=x+h;
yx=y(x);
printf("%3.1f %9.6f %9.6f %9.6f \n"
,x,ye,ym,yr);
}
}
MATLAB版本
f.m
function f = f(x,y)
f = ((y)-2*(x)/(y));
end
y.m
function y = y( x )
y = sqrt(1+2*(x));
end
主程序（记得把上面两个文件放在同级目录，我真是操碎了心）
close all;clc;clear;
a=0;b=0.9;h=0.1;y0=1;
x=a;
ye=y0;
ym=y0;
yr=y0;
while(x<=b)
ye=ye+h*f(x,ye);
yp=ym+h*f(x,ym);
ym=ym+h/2*(f(x,ym)+f(x+h,yp));
k1=f(x,yr);
k2=f(x+h/2,yr+h/2*k1);
k3=f(x+h/2,yr+h/2*k2);
k4=f(x+h,yr+h*k3);
yr=yr+h/6*(k1+2*k2+2*k3+k4);
x=x+h;
yx=y(x);
fprintf(' %.6f %.6f %.6f\n',ye,ym,yr)
end
sprintf('三列分别为欧拉方法、改进的欧拉方法、四阶龙格库塔方法\n每行对应的x从0.1开始以h=0.1跃进到1,共10行')
记录两个小坑
牛顿插值法的MATLAB版本里有一个下标是j-i+2，达到的目的是下标非负取整，具体逻辑要想清楚
要用getchar把input的变为数字，不然是字符会有问题，使用str2double也行不过一般是getchar
最后的一些想法
估计可爱的学弟学妹copy完代码就直接关了，有能力的也不稀罕这些不是很专业的代码。
c语言版本的是从一本书上参考的（忘了叫啥名字了，但挺有年代感的是零几年的，在学校图书馆里吃灰），MATLAB版本是我照着c语言版本自己“翻译”的，水平明显比不上，但都可以跑的通可以正确运行。
再多说点，好像是在一个问答里看到的，说，数学系的同学要具备代码能力，这种能力不是说要像计科的同学一样，考虑复杂度啊代码可读性啊换行啊什么特别专业的，但是处理数据完成你想要完成的能力还是必须要有的，希望数学系同志们在未来的日子里共同持续进化
共勉
