> | restart; SecDeri := proc(x,y,h,n) local i,j,A,b; with(linalg): A:=array(1..n,1..n); b:=array(1..n); # Reset A to 0 for i from 1 to n do for j from 1 to n do A[i,j]:=0; od; od; A[1,1]:=1; b[1]:=0; A[n,n]:=1; b[n]:=0; for i from 2 to n-1 do A[i,i-1]:=h[i-1]/6; A[i,i]:=(h[i-1]+h[i])/3; A[i,i+1]:=h[i]/6; b[i]:=(y[i+1]-y[i])/h[i]-(y[i]-y[i-1])/h[i-1]; od; RETURN(linsolve(A,b)); end: CubicSpline := proc() local fp,filename, N1,X1,Y1, N2,X2,Y2, F1,F2,P1,P2, i,j, H1,M1, H2,M2, a,b,c,d, Udp,Ldp,U,L; with(plots): printf("***** Author: Mengnien Wu *************************************************\n"); printf("* Draw 2 cubic splines for N1 and N2 data points respectively. *\n"); printf("* Data points are supposed obtained from a graph that can be decomposed *\n"); printf("* into 2 parts, upper and lower, graphs of some continuous functions, and *\n"); printf("* N1 points from the upper part and N2 from the lower one are picked. *\n"); printf("* A file which contains *\n"); printf("* N1, X1[1],Y1[1] ... X1[N1],Y1[N1], *\n"); printf("* N2, X2[1],Y2[1] ... X2[N2],Y2[N2] *\n"); printf("* will be read. *\n"); printf("***************************************************************************\n"); printf("\nEnter your filename:\n"); filename := scanf("%a")[1]; fp:=fopen(filename,READ); N1:=fscanf(fp,"%d")[1]; X1:=array(1..N1); Y1:=array(1..N1); for i from 1 to N1 do X1[i]:=fscanf(fp,"%f")[1]; Y1[i]:=fscanf(fp,"%f")[1]; od; N2:=fscanf(fp,"%d")[1]; X2:=array(1..N2); Y2:=array(1..N2); for i from 1 to N2 do X2[i]:=fscanf(fp,"%f")[1]; Y2[i]:=fscanf(fp,"%f")[1]; od; fclose(fp); F1:=array(1..N1-1); #Function P1:=array(1..N1-1); #Picture H1:=array(1..N1-1); for i from 1 to N1-1 do H1[i]:=X1[i+1]-X1[i]; od; M1:=SecDeri(X1,Y1,H1,N1); for j from 1 to N1-1 do a :=Y1[j]; b :=(Y1[j+1]-Y1[j])/H1[j]-(2*M1[j]+M1[j+1])*H1[j]/6; c :=M1[j]/2; d :=(M1[j+1]-M1[j])/(6*H1[j]); F1[j]:=a+(x-X1[j])*( b+(x-X1[j])*( c+(x-X1[j])*d ) ); #Horner's P1[j]:=plot(F1[j],x=X1[j]..X1[j+1],color=blue,scaling=constrained): od; F2:=array(1..N2-1); #Function P2:=array(1..N2-1); #Picture H2:=array(1..N2-1); for i from 1 to N2-1 do H2[i]:=X2[i+1]-X2[i]; od; M2:=SecDeri(X2,Y2,H2,N2); for j from 1 to N2-1 do a := Y2[j]; b := (Y2[j+1]-Y2[j])/H2[j] - (2*M2[j]+M2[j+1])*H2[j]/6; c := M2[j]/2; d := (M2[j+1]-M2[j])/(6*H2[j]); F2[j]:= a+(x-X2[j])*( b+(x-X2[j])*( c+(x-X2[j])*d ) ); #Horner's P2[j]:= plot(F2[j],x=X2[j]..X2[j+1],color=blue,scaling=constrained): od; Udp:=plot([ [X1[n],Y1[n]] $n=1..N1 ],style=point,symbol=circle): Ldp:=plot([ [X2[n],Y2[n]] $n=1..N2 ],style=point,symbol=circle): display(Udp,Ldp,[P1[n] $n=1..N1-1],[P2[n] $n=1..N2-1]); end: |
> | CubicSpline(); |
Warning, the name changecoords has been redefined
***** Author: Mengnien Wu *************************************************
* Draw 2 cubic splines for N1 and N2 data points respectively. *
* Data points are supposed obtained from a graph that can be decomposed *
* into 2 parts, upper and lower, graphs of some continuous functions, and *
* N1 points from the upper part and N2 from the lower one are picked. *
* A file which contains *
* N1, X1[1],Y1[1] ... X1[N1],Y1[N1], *
* N2, X2[1],Y2[1] ... X2[N2],Y2[N2] *
* will be read. *
***************************************************************************
Enter your filename:
> | "DataPoints.txt" |
Warning, the protected names norm and trace have been redefined and unprotected
> |