FORtran多项式插值.

FORTRAN 90 程序编程规范

FORTRAN 90 程序编程规范 Fortran 90 编程规范，使程序代码高度组织化，更加易读、易懂、易于维护，程序更加高效。使编出的程序更易懂、易于维护。 1 语言选择 数值预报创新系统软件开发应避免使用Fortran77 的某些过时特征以Fortran 90不一致的特征。选择Fortran 90 作为开发语言，并采用Fortran 90 的新功能，如动态内存的分配(dynamic memory allocation)、递归（recursion ）, 模块(modules)、POINTER 、长变量名、自由格式等。 Fortran 77其中某些只是一些冗余的功能，这些功能已经过时，另外，还有一些在Fortran90 中被证明是不好的用法，建议不要使用。 2 Fortran 90 的新特性 2.1.1 建议使用的Fortran 90 新特性 建议使用Fortran 90 提供的模块（module ），并用Use ONLY 指定module 中哪些变量或派生类型定义可用于调用程序。 尽量使用数组下标三元组，这样可优化并减少所需的代码行数。为提高可读性，要在括号内表明数组的维数，例如： 1dArrayA(:) = 1dArrayB(:) + 1dArrayC(:) 2dArray(: , :) = scalar * Another2dArray(: , :) 当访问数组的子集时，例如在有限差分等式中，可以通过使用下标三元组实现。例如：2dArray(: , 2:len2) = scalar *（ & Another2dArray(:, 1:len2 -1) & - Another2dArray(:, 2:len2) & ) 对程序单元（program units ）命名，并使用End program ，End subroutine ，End interface ，End module 等结构再次指定“program unit ”的名称。 在逻辑表达式中使用>、 >=、 ==、 <、 <=、 /=，它们分别代 替.gt.、.ge.、.eq.、.lt.、.le.、.ne. 。新的表示方法更接近标准的数学符号 在变量定义中始终使用“：：”；始终用“DIMENSION ”定义数组形状；始终用（len=）的语法格式声明字符变量的长度。

三次样条插值代码

2 三次样条插值程序 三次样条插值利用方案二（求解固支样条或压紧样条） 按照要求要起点和终点的一阶导数值已知， 可得关于01,,.....,n M M M 的严格对角占优势的三对角方程组 然后利用三对角法（追赶法）解此线性方程组。 （1）编写M 文件，并保存文件名scfit.m % x,y 分别为n 个节点的横坐标和纵坐标值组成的向量 % dx0和dxn 分别为S 的导数在x0和xn 处的值，即m 0和m n n=length(x)-1; h=diff(x); d=diff(y)./h; a=h(2:n-1); b=2*(h(1:n-1)+h(2:n)); c=h(2:n); u=6*diff(d); b(1)=b(1)-h(1)/2; u(1)=u(1)-3*(d(1)-dx0); b(n-1)=b(n-1)-h(n)/2; u(n-1)=u(n-1)-3*(dxn-d(n)); %追赶法部分 for k=2:n-1 temp=a(k-1)/b(k-1); b(k)=b(k)-temp*c(k-1); u(k)=u(k)-temp*u(k-1); end m(n)=u(n-1)/b(n-1); for k=n-2:-1:1 m(k+1)=(u(k)-c(k)*m(k+2))/b(k); end %求S K1，S K2，S K3，S K4 m(1)=3*(d(1)-dx0)/h(1)-m(2)/2; m(n+1)=3*(dxn-d(n))/h(n)-m(n)/2; for k=0:n-1 00 ()S x m '=()n n S x m '=0011111111212212n n n n n n M d M d M d M d μλμλ----??????????????????????=??????????????????????????

fortran程序案例题汇编(14道)

1．Fibonacci数列定义如下： F1=1 F2=1 F n=F n-1+F n-2 (n>2) 求Fibonacci数列的前30项。 integer f(30),i f(1)=1 f(2)=2 do i=3,30 f(i)=f(i-1)+f(i-2) enddo print*,f end 2．输入10个学生的总分，求每个学生的名次integer s(10),a(10),i,j do i=1,10 read*,s(i) enddo do i=1,10 a(i)=1 do j=1,10 if(s(i)

3．给定一组数，按照从小到大的顺序输出。integer a(10) integer p do i=1,10 read *,a(i) enddo do j=1,9 p=j do i=j+1,10 if (a(i)

4．输入若干名学生的学号和三门课程（语数英）的成绩，要求从键盘输入一个学生的学号，能打印出该学生的三门功课成绩和总分。 character*6,dimension(:),allocatable::xue integer,dimension(:,:),allocatable::g integer,dimension(:),allocatable::zong integer i,j,n character*6,xh print *,"请输入学生的个数" read *,n allocate(xue(n)) allocate(g(n,3)) allocate(zong(n)) do i=1,n read *,xue(i),(g(i,j),j=1,3) enddo do i=1,n zong(i)=0 do j=1,3 zong(i)=zong(i)+g(i,j) enddo enddo print *,"请输入你要打印的学生的学号" read *,xh do i=1,n if(xue(i)==xh)then print *,(g(i,j),j=1,3),zong(i) exit endif enddo end

三次样条插值函数

沈阳航空航天大学 数学软件课程设计 （设计程序） 题目三次样条插值函数 班级 / 学号 学生姓名 指导教师

沈阳航空航天大学 课程设计任务书 课程名称数学软件课程设计 院（系）理学院专业信息与计算科学 班级学号姓名 课程设计题目三次样条插值函数 课程设计时间: 2010 年12月20日至2010 年12月31日 课程设计的内容及要求： 1.三次样条插值函数 给出函数在互异点处的值分别为。 (1)掌握求三次样条插值函数的基本原理； (2)编写程序求在第一边界条件下函数的三次样条插值函数； (3)在区间上取n=10，20，分别用等距节点对函数 作三次样条插值函数，利用(1)的结果画出插值函数的图形，并在该图形界面中同时画出的图形。 [要求] 1.学习态度要认真，要积极参与课程设计，锻炼独立思考能力； 2.严格遵守上机时间安排； 3.按照MATLAB编程训练的任务要求来编写程序； 4.根据任务书来完成课程设计论文； 5.报告书写格式要求按照沈阳航空航天大学“课程设计报告撰写规范”； 6.报告上交时间：课程设计结束时上交报告；

7.严谨抄袭行为。 指导教师年月日负责教师年月日学生签字年月日

沈阳航空航天大学 课程设计成绩评定单 课程名称数学软件课程设计 院（系）理学院专业信息与计算科学课程设计题目三次样条插值函数 学号姓名 指导教师评语： 课程设计成绩 指导教师签字 年月日

目录 一正文 (1) 1问题分析 (1) 1.1 题目 (1) 1.2 分析 (1) 2 研究方法原理 (1) 2.1 求三次样条插值多项式，算法组织 (1) 3 算例结果 (3) 二总结 (7) 参考文献 (8) 附录 (9) 源程序: (9) 程序1 (9) 程序2 (10) 程序3 (12) 程序 4 (12)

三次样条插值方法的应用

CENTRAL SOUTH UNIVERSITY 数值分析实验报告

三次样条插值方法的应用 一、问题背景 分段低次插值函数往往具有很好的收敛性，计算过程简单，稳定性好，并且易于在在电子计算机上实现，但其光滑性较差，对于像高速飞机的机翼形线船体放样等型值线往往要求具有二阶光滑度，即有二阶连续导数，早期工程师制图时，把富有弹性的细长木条（即所谓的样条）用压铁固定在样点上，在其他地方让他自由弯曲，然后沿木条画下曲线，称为样条曲线。样条曲线实际上是由分段三次曲线并接而成，在连接点即样点上要求二阶导数连续，从数学上加以概括就得到数学样条这一概念。下面我们讨论最常用的三次样条函数及其应用。 二、数学模型 样条函数可以给出光滑的插值曲线（面），因此在数值逼近、常微分方程和偏微分方程的数值解及科学和工程的计算中起着重要的作用。 设区间[]b ,a 上给定有关划分b x x n =<<<= 10x a ，S 为[]b ,a 上满足下面条件的函数。 ● )(b a C S ,2∈； ● S 在每个子区间[]1,+i i x x 上是三次多项式。 则称S 为关于划分的三次样条函数。常用的三次样条函数的边界条件有三种类型： ● Ⅰ型 ()()n n n f x S f x S ''0'',==。 ● Ⅱ型 ()()n n n f x S f x S ''''0'''',==，其特殊情况为()()0''''==n n x S x S 。 ● Ⅲ型 ()() 3,2,1,0,0==j x S x S n j j ，此条件称为周期样条函数。 鉴于Ⅱ型三次样条插值函数在实际应用中的重要地位，在此主要对它进行详细介绍。 三、算法及流程 按照传统的编程方法，可将公式直接转换为MATLAB 可是别的语言即可；另一种是运用矩阵运算，发挥MATLAB 在矩阵运算上的优势。两种方法都可以方便地得到结果。方法二更直观，但计算系数时要特别注意。这里计算的是方法一的程序，采用的是Ⅱ型边界条件，取名为spline2.m 。 Matlab 代码如下： function s=spline2(x0,y0,y21,y2n,x) %s=spline2(x0,y0,y21,y2n,x) %x0,y0 are existed points,x are insert points,y21,y2n are the second

fortran程序实例

1)实例3—求多个半径下的圆周长 ! z3.f90 --Fortran95 ! FUNCTIONS: ! z3 - Entry point of console application. !************************************************************************* ! PROGRAM: z3 ! PURPOSE: Entry point for the console application. !************************************************************************ program z3 ! 求多个半径下的圆周长 ! 主程序 ! PROGRAM Z3 PRINT *, 'R=',1.2,'C=',C(1.2) PRINT *, 'R=',3.4,'C=',C(3.4) PRINT *, 'R=',15.6,'C=',C(15.6) PRINT *, 'R=',567.3,'C=',C(567.3) END program z3 !子程序 FUNCTION C(R) PI=3.1415926 C=2*PI*R RETURN ! Body of z3 end

2)实例4—键盘与显示器输入/输出 a)Fortran 基本操作 b)程序指令 ! ZXZ_I_O.f90 ! FUNCTIONS: ! ZXZ_I_O - Entry point of console application. ! PROGRAM: ZXZ_I_O ! PURPOSE: Entry point for the console application. !***************输入、输出样式种种************************** program ZXZ_I_O implicit none !变量声明的位置 INTEGER(2) i; INTEGER(4) j; INTEGER(4) m; REAL n INTEGER A,B ! Variables PRINT*,'输入整数A'; READ*, A PRINT*,'输入整数B'; READ*, B B=A+B PRINT*,'B=A+B=',B WRITE(*,*) 'A*B=',A*B PRINT* ,'以上为计算机的计算结果，注意B的值'

数值分析作业-三次样条插值

数值计算方法作业 实验4.3 三次样条差值函数 实验目的： 掌握三次样条插值函数的三弯矩方法。 实验函数: dt e x f x t ? ∞ -- = 2 221)(π 实验内容： （1） 编程实现求三次样条插值函数的算法，分别考虑不同的边界条件； （2） 计算各插值节点的弯矩值； （3） 在同一坐标系中绘制函数f(x)，插值多项式，三次样条插值多项式的曲线 比较插值结果。 实验4.5 三次样条差值函数的收敛性 实验目的： 多项式插值不一定是收敛的，即插值的节点多，效果不一定好。对三次样条插值函数如何呢？理论上证明三次样条插值函数的收敛性是比较困难的，通过本实验可以证明这一理论结果。 实验内容： 按照一定的规则分别选择等距或非等距的插值节点，并不断增加插值节点的个数。 实验要求： （1） 随着节点个数的增加，比较被逼近函数和三样条插值函数的误差变化情 况，分析所得结果并与拉格朗日插值多项式比较； （2） 三次样条插值函数的思想最早产生于工业部门。作为工业应用的例子，考

虑如下例子：某汽车制造商根据三次样条插值函数设计车门曲线，其中一 算法描述： 拉格朗日插值： 错误!未找到引用源。 其中错误!未找到引用源。是拉格朗日基函数，其表达式为：() ∏ ≠=--=n i j j j i j i x x x x x l 0) ()( 牛顿插值： ) )...()(](,...,,[.... ))(0](,,[)0](,[)()(1102101210100----++--+-+=n n n x x x x x x x x x x f x x x x x x x f x x x x f x f x N 其中????? ?? ?? ?????? --=--= --= -)/(]),...,[],...,[(]...,[..],[],[],,[)()(],[01102110x x x x x f x x x f x x x f x x x x f x x f x x x f x x x f x f x x f n n n n i k j i k j k j i j i j i j i 三样条插值： 所谓三次样条插值多项式Sn(x)是一种分段函数，它在节点Xi(a

几个fortran程序

计算圆周率 REAL R,R1,R2,PI ISEED=RTC() N0=0 N=300000 DO I=1,N R1=RAN(ISEED) R2=RAN(ISEED) R=SQRT(R1*R1+R2*R2) IF(R<1.0)N0=N0+1 END DO PI=4.0*N0/N WRITE(*,*)PI END 一）蒙特卡洛计算生日问题 假设有N个人在一起，各自的生日为365天之一，根据概率理论，与很多人的直觉相反，只需23个人便有大于50％的几率人群中至少有2个人生日相同。 INTEGER M(1:10000), NUMBER1(0:364), NUMBER2 REAL X,Y ISEED=RTC() DO J=1, 10000 NUMBER1=0 X=RAN(ISEED) NUMBER1(0)=INT(365*X+1) JJJ=1 DO I=1,364 Y=RAN(ISEED) NUMBER2=INT(365*Y+1) ETR=COUNT(NUMBER1.EQ.NUMBER2) IF (ETR= =1) THEN EXIT ELSE JJJ=JJJ+1 M(J)=JJJ NUMBER1(I)=NUMBER2 END IF END DO END DO DO I=1,10000 IF(M(I).LE.23) SUM=SUM+1 END DO PRINT *,SUM/10000 END 二）MONTE CARLO SIMULATION OF ONE DIMENSIONAL DIFFUSION 蒙特卡罗计算一维扩散问题

INTEGER X,XX(1:1000,1:1000) REAL XXM(1:1000) ! X:INSTANTANEOUS POSITION OF ATOM ! XX(J,I)：X*X ,J:第几天实验,I:第几步跳跃 ! XXM(I): THE MEAN OF XX WRITE(*,*) "实验天数JMAX，实验次数IMAX" READ(*,*) JMAX,IMAX ISEED=RTC() DO J=1,JMAX !第几天实验 X=0 !!! DO I=1,IMAX !第几步跳跃 RN=RAN(ISEED) IF(RN<0.5)THEN X=X+1 ELSE X=X-1 END IF XX(J,I)=X*X END DO END DO OPEN(1,FILE="C:\DIF1.DAT") DO I=1,IMAX XXM=0.0 XXM(I)=1.0*SUM(XX(1:JMAX,I))/JMAX !! WRITE(1,*) I, XXM(I) END DO CLOSE(1) END 三维的！ 三）通过该程序了解FORTRAN语言如何画图（通过像素画图）USE MSFLIB INTEGER XR，YR !在的区域中画一个圆 PARAMETER XR=400,YR=400 INTEGER R, S(1:XR,1:YR) X0=XR/2 ! 圆心位置X0，YO Y0=YR/2 R=MIN(X0-10,Y0-10) !圆半径 S=0 !像素的初始状态（颜色） DO I=1,XR DO J=1,YR IF((I-X0)**2+(J-Y0)**2<=R**2)S(I,J)=10 IER=SETCOLOR(S(I,J)) IER=SETPIXEL(I,J) END DO END DO END

fortran程序30个

程序1 PROGRAM chaper1272 INTEGER :: A(10)=(/5,7,4,8,12,2,10,3,9,11/) INTEGER :: sum=0,ave !打开一数据文件,设置一个有格式直接存取文件,将10个数分2个记录写入文件。 !数据文件生成2个记录,每个记录行长度相同,记录长度为25 OPEN(1,FILE='input21.dat',FORM='FORMATTED',ACCESS='D IRECT',RECL=25) WRITE(1,"(5I5)",REC=1)(A(I)+10,I=1,5) !按格式说明将头5个数写入第1个记录 WRITE(1,"(5I5)",REC=2)(A(I)+10,I=6,10)!按格式说明将后5个数写入第2个记录 READ(1,"(5I5)",REC=2)(A(I),I=6,10) !按格式说明从第2个记录中读取后5个数 READ(1,"(5I5)",REC=1)(A(I),I=1,5) !按格式说明从第1个记录中读取头5个数 DO I=1,10 sum=sum+A(I) ENDDO ave=sum/10 !打开一个最大记录长度为22的有格式顺序存取文件

OPEN(2,FILE='input22.dat',FORM='FORMATTED',ACCESS='D IRECT',RECL=22) WRITE(2,"('10个数之和为：',I5)",REC=1) sum !输出1记录行,记录长度为22 WRITE(2,"('10个数平均值为：',I5)",REC=2)ave !输出1记录行,记录长度为22 WRITE(2,"(A)",REC=3) '程序运行正常结束。' !输出一个记录行,记录长度为22 END 程序2 PROGRAM average REAL sum,ave INTEGER n OPEN(1,file='score.dat') PRINT*,'正在统计平均成绩,请等待。' sum=0.0;n=0 DO READ(1,*,END=100) S sum=sum+s n=n+1 ENDDO 100 ave=sum/n

(精选)三次样条插值的MATLAB实现

MATLAB 程序设计期中考查 在许多问题中，通常根据实验、观测或经验得到的函数表或离散点上的信息，去研究分析函数的有关特性。其中插值法是一种最基本的方法，以下给出最基本的插值问题——三次样条插值的基本提法： 对插值区间[]b a ,进行划分：b x x x a n ≤

关于三次样条插值函数的学习报告(研究生)资料

学习报告—— 三次样条函数插值问题的讨论 班级：数学二班 学号：152111033 姓名：刘楠楠

样条函数： 由一些按照某种光滑条件分段拼接起来的多项式组成的函数；最常用的样条函数为三次样条函数，即由三次多项式组成，满足处处有二阶连续导数。 一、三次样条函数的定义： 对插值区间[,]a b 进行划分，设节点011n n a x x x x b -=<< <<=，若 函数2()[,]s x c a b ∈在每个小区间1[,]i i x x +上是三次多项式，则称其为三次样条函数。如果同时满足()()i i s x f x = (0,1,2)i n =，则称()s x 为()f x 在 [,]a b 上的三次样条函数。 二、三次样条函数的确定： 由定义可设：101212 1(),[,] (),[,]()(),[,] n n n s x x x x s x x x x s x s x x x x -∈??∈?=???∈?其中()k s x 为1[,]k k x x -上的三次 多项式，且满足11(),()k k k k k k s x y s x y --== (1,2,,k n = 由2()[,]s x C a b ∈可得：''''''()(),()(),k k k k s x s x s x s x -+-+== 有''1()(),k k k k s x s x -++= ''''1()(),(1 ,2,,1)k k k k s x s x k n -+ +==-， 已知每个()k s x 均为三次多项式，有四个待定系数，所以共有4n 个待定系数，需要4n 个方程才能求解。前面已经得到22(1)42n n n +-=-个方程，因此要唯一确定三次插值函数，还要附加2个条件，一般上，实际问题通常对样条函数在端点处的状态有要求，即所谓的边界条件。 1、第一类边界条件：给定函数在端点处的一阶导数，即 ''''00(),()n n s x f s x f == 2、第二类边界条件：给定函数在端点处的二阶导数，即

Fortran95程序设计课后习题答案(word版方便)

第四章 1.program main implicit none write(*,*) "Have a good time." write(*,*) "That's not bad." write(*,*) '"Mary" isn''t my name.' end program 2.program main real, parameter :: PI=3 implicit none.14159 real radius write(*,*) "请输入半径长" read(*,*) radius write(*,"(' 面积='f8. 3)") radius*radius*PI end program 3.program main implicit none real grades write(*,*) "请输入成绩" read(*,*) grades write(*,"(' 调整后成绩为'f8.3)") SQRT(grades)*10.0 end program 4.integer a,b real ra,rb a=2 b=3 ra=2.0 rb=3.0 write(*,*) b/a ! 输出1, 因为使用整数计算, 小数部分会无条件舍去write(*,*) rb/ra ! 输出1.5 5.p rogram main implicit none type distance real meter, inch, cm end type type(distance) :: d write(*,*) "请输入长度:" read(*,*) d%meter d%cm = d%meter*100 d%inch = d%cm/2.54 write(*,"(f8.3'米='f8.3'厘米='f8.3'英寸')") d%meter, d%cm, d%inch end program 第五章 1.program main implicit none integer money real tax write(*,*) "请输入月收入" read(*,*) money if ( money<1000 ) then tax = 0.03 else if ( money<5000) then tax = 0.1 else tax = 0.15 end if write(*,"(' 税金为'I8)") nint(money*tax) end program 2.program main implicit none integer day character(len=20) :: tv write(*,*) "请输入星期几" read(*,*) day select case(day) case(1,4) tv = "新闻" case(2,5) tv = "电视剧" case(3,6) tv = "卡通" case(7) tv = "电影" case default write(*,*) "错误的输入" stop end select write(*,*) tv end program 3.program main implicit none integer age, money real tax write(*,*) "请输入年龄" read(*,*) age write(*,*) "请输入月收入" read(*,*) money if ( age<50 ) then if ( money<1000 ) then tax = 0.03 else if ( money<5000 )then tax = 0.10 else tax = 0.15 end if else if ( money<1000 ) then tax = 0.5 else if ( money<5000 )then tax = 0.7 else tax = 0.10 end if end if write(*,"(' 税金为'I8)") nint(money*tax) end program 4.program main implicit none integer year, days logical mod_4, mod_100, mod_400 write(*,*) "请输入年份" read(*,*) year mod_4 = ( MOD(year,4) == 0 ) mod_100 = ( MOD(year,100) == 0 ) mod_400 = ( MOD(year,400) == 0 ) if ( (mod_4 .NEQV. mod_100) .or. mod_400 ) then days = 366 else days = 365 end if write(*,"('这一年有'I3'天')") days stop end program 第六章 1.program main implicit none integer i do i=1,5 write(*,*) "Fortran" end do stop end program 2.program main implicit none integer i,sum sum = 0 do i=1,99,2 sum = sum+i end do write(*,*) sum stop end program 3.program main implicit none integer, parameter :: answer = 45 integer, parameter :: max = 5 integer weight, i do i=1,max write(*,*) "请输入体重" read(*,*) weight if ( weight==answer ) exit end do if ( i<=max ) then write(*,*) "猜对了" else write(*,*) "猜错了" end if stop end program 4.program main implicit none integer, parameter :: max=10 integer i real item real ans ans = 1.0 item = 1.0 do i=2,max item = item/real(i) ans = ans+item

FORTRAN经典入门程序20例

对于FORTRAN的初学者。这些例子可作为小练习。 1.例题：计算工资问题。每小时工资为RATE,如果工作超过40小时，加班呢部分工资是正常时间工资的1.5倍。 C Payroll with overtime program payroll real rate, hours, pay read (*,*) rate, hours if (hours>40) then pay=40*rate+(hours-40)*1.5*rate else pay=hours*rate END IF print *,"rate=", rate print *, "hours=", hours print *,"pay=",pay end 2.学生成绩问题。大于80为A级。大于60小于80为B级。小于60为C级。IF的嵌套。注意空格可以看清楚else if ,end if,print的内容. PROGRAM GRADESTUDENT READ *,GRADE IF (GRADE .GE. 80) THEN PRINT *,GRADE,"GRADE=>A" ELSE IF (GRADE.LT.60) THEN PRINT*,GRADE,"GRADE=>C" ELSE PRINT*,GRADE,"GRADE=>B" END IF END IF END 3.三个数按从小到大排序。 PROGRAM MAXMIN REAL A,B,C,T READ *,A,B,C IF (A.GT.B) THEN T=A A=B B=T ELSE END IF IF (B.GT.C) THEN T=B B=C

PROGRAM EQUATION READ (*,*) X IF (X .GE. 0.0) Y=3*X+6 IF (X .LT. 0.0) Y=-X**2+2*X-8 PRINT *,"X=",X, "Y=",Y END 6.CONTINUE语句。为了使循环的起止范围清晰，使用终端语句 例：求5！ program ex1 FACT=1 do 10, N=1,5,1 ##DO 语句标号，循环不变量=初值，终值，步长FACT=FACT*N 10continue print*,FACT end 求1！+2！+3！+···+N！ READ *,N FACT=1

MATLAB三次样条插值之三弯矩法

ＭATLＡＢ三次样条插值之三弯矩法 首先说这个程序并不完善，为了实现通用（１，２，…,ｎ）格式解题，以及为调用追赶法程序，没有针对节点数在三个以下的情况进行分类讨论。希望能有朋友给出更好的方法。 首先，通过函数saｎwanｊ得到方程的系数矩阵，即追赶法方程的四个向量参数,接下来调 用追赶法(在inteｒｓaｎwj函数中），得到三次样条分段函数系数因子,然后进行多项式合并 得到分段函数的解析式，程序最后部分通过判断输入值的区间自动选择对应的分段函数并计算 改点的值。附:追赶法程序cｈase ％%%％％%%%%%%%％% fｕncｔiｏｎ [nｅwv,w,newu，neｗd]=sanwj（ｘ,y,x０,y0,y1ａ,y1b）?％三弯矩样 条插值?％将插值点分两次输入，x0y０单独输入?% 边值条件a的二阶导数 y1a 和ｂ 的二阶导数ｙ1b，这里建议将y１ａ和y1b换成y2a和ｙ2b,以便于和三转角代码相区别 ?n=length(x);m=leｎgth(y); ｉf m~=n?errｏｒ(＇x or y 输入有误,再来'）; end?v=ｏneｓ(ｎ-１，1);u=ｏnes(n－１,1）;d=zeros（n－１,１）;?w=２*ｏ nes(n+１)；?h0=x（１）-x0;?h=zerｏs(n-１,1); for k＝1:n-１?ｈ(k)＝x（k+1)－x（k);?enｄ ｖ(1)＝h0/(ｈ0+h(1）); u(1)＝1－v（1）; ｄ（１)＝6*((y（2）-y(1）)/h（1)-(y(1)-y０)/ｈ0）/(h0+ｈ(１）);?% ｆoｒ k＝2:ｎ-1?ｖ(ｋ）＝h(ｋ-１)/(ｈ(k-１)+h(ｋ)）;?u(k)=1－v(k);?d（k)＝ 6*((y(ｋ+１)-ｙ（ｋ)）/h（k）-(ｙ(k）-y(ｋ-1）)/ｈ(k－1))／(h（k-1)+h(k)); ｅnd newv=[v;1］;?newu=［1;u]； d0=6＊（(y(1)-y0)/h0-y1a)/h0; d(n)=6＊（y1ｂ-(ｙ(n）－ｙ(n-1))／ｈ(n-1))／h(ｎ-1）; newd=［d0;d]; %%%%%%％％%％％％ ｆunctｉon iｎｔersａnｗj（x，y,ｘ0,y0,y1a，y1b) %三弯矩样条插值?%第一部分?n=leｎgtｈ(ｘ）;ｍ＝length(y); iｆ m~＝n?eｒror（'ｘｏｒy 输入有误，再来'); enｄ?%重新定义h?h=zｅｒｏs(n,1）; h(1)=x（1）-ｘ０; fｏr k=2:n h（k)=ｘ(ｋ)－x(k-1)；?end ％sptep1调用三弯矩函数?[a,b,c,d]=sanwj(x，ｙ,ｘ0,ｙ0,y１ａ,y1b);

优秀的fortran程序编程规范

Programming Guidelines for PARAMESH Software Development(NOTE: This document is heavily based upon the Introduction This document describes the programming guidelines to be used by software developers wishing to contribute software to the PARAMESH, parallel, adaptive mesh refinement software. We welcome people to contribute software and/or bug fixes to the PARAMESH AMR software. Software to be added to PARAMESH can come in 2 forms: ?Improvements to the basic PARAMESH kernal software found in the mpi_source, source and hearders directories. ?Software the addes additional functionality to PARAMESH. This type of software should be added as separate entities within the utilities directory. Complete applications should not be added as part of PARAMESH. PARAMESH is only meant to be a tool which supports parallel adaptive mesh applications and any software which supports this goal will be considered for acceptance into PARAMESH. For instance, a solver for the poisson equation that works with PARAMESH would be acceptable, but an application that solves the equation of gas dynamics would not. The PARAMESH software is slowly being evolved to be consistent with this document. Any new software which is contributed should follow these guidlines. If not, it will be rejected. This document deals mainly with Fortran 90, since most new PARAMESH software will probably be written in that language. [Throughout this document, the term "Fortran" should be understood to mean Fortran 90.] Since we expect C and C++ also to be used, a separate document dealing with them will be developed in the future. In the meantime, this document can serve as a general guideline for developing code to be used with PARAMESH in those programming languages. The guidlines in this document should be adhered to by ANY software which will be released as part of the PARAMESH package of source code. This includes software 'utilities' (stored in the paramesh/utilities directory) which add functionality to PARAMESH for different algorithms. It also should be applied to any new code developed and added to the main source code for PARAMESH in the paramesh/source, paramesh/mpi_source, or paramesh/headers directories.

平面四边形四节点等参单元Fortran源程序

C ************************************************ C * FINITE ELEMENT PROGRAM * C * FOR Two DIMENSIONAL ELASticity PROBLEM * C * WITH 4 NODE * C ************************************************ PROGRAM ELASTICITY character*32 dat,cch DIMENSION SK(80000),COOR(2,300),AE(4,11),MEL(5,200), & WG(4),JR(2,300),MA(600),R(600),iew(30),STRE(3,200) COMMON /CMN1/ NP,NE,NM,NR COMMON /CMN2/ N,MX,NH COMMON /CMN3/ RF(8),SKE(8,8),NN(8) WRITE(*,*)'PLEASE ENTER INPUT FILE NAME' READ(*,'(A)')DAT OPEN(4,FILE=dat,STATUS='OLD') OPEN(7,FILE='OUT',STATUS='UNKNOWN') READ(4,*)NP,NE,NM,NR WRITE(7,'(A,I6)')'NUMBER OF NODE---------------------NP=',np WRITE(7,'(A,I6)')'NUMBER OF ELEMENT------------------NE=',ne WRITE(7,'(A,I6)')'NUMBER OF MATERIAL-----------------NM=',nm WRITE(7,'(A,I6)')'NUMBER OF surporting---------------NC=',Nr CALL INPUT (JR,COOR,AE,MEL)