Форум сайта python.su
from numpy import * import pylab from matplotlib import mlab from scipy.linalg import solve from scipy import array import matplotlib.pyplot as plt from pylab import * f1=array([ [4.0,102.56], [4.2,104.18], [4.5,100.11]]) def sweep(f,d2fa,d2fb): k=2 A = zeros(len(f)) B = zeros(len(f)) C = zeros(len(f)) F = zeros(len(f)) res = zeros(len(f)) alpha = zeros(len(f)) beta = zeros(len(f)) while (k<len(f)): A[k]=h(f,k-1) B[k]=h(f,k) C[k]=2*(h(f,k-1)+h(f,k)) F[k]=3*(func(f,k)-func(f,k-1)) k+=1 alpha[2]=-B[2]/C[2] beta[2]=F[2]/C[2] k=2 while (k<len(f)-1): alpha[k+1]=-B[k]/(A[k]*alpha[k]+C[k]) beta[k+1]=(F[k]-A[k]*beta[k])/(A[k]*alpha[k]+C[k]) k+=1 k=len(f)-2 res[0]=d2fa res[len(f)-1]=d2fb while(k>0): res[k]=alpha[k+1]*res[k+1]+beta[k+1] k-=1 return res def h(f,k): return f[k,0]-f[k-1,0] def func(f,k): return (f[k,1]-f[k-1,1])/h(f, k) def sweep2(f,d2fa,d2fb): k=2 sigma = zeros(10) lamda = zeros(10) C=zeros(10) C[0]=d2fa C[9]=d2fb sigma[0]=-h(f, 2)/(2*h(f,1)+2*h(f,2)) lamda[0]=3*(func(f,2)-func(f,1))/(2*(h(f,1)+h(f,2))) while k<len(f): sigma[k-1]=h(f,k)/(2*h(f,k-1)+2*h(f,k)+h(f,k-1)*sigma[k-2]) lamda[k-1]=(3*(func(f,k)-func(f,k-1))-h(f,k-1)*lamda[k-2])/(2*h(f,k-1)+2*h(f,k)+h(f,k-1)*sigma[k-2]) k+=1 k=9 while k>1: C[k-1]=sigma[k-1]*C[k]+lamda[k-1] k-=1 return C def ispline(f,x,C): if x<f[0,0]: print "viberite bol'she x" elif x>f[len(f)-1,0]: print "viberite men'she x" else: k=0 count=0 A = f[:,1] B = zeros(len(f)) D = zeros(len(f)) while (k<len(f)): if (x-f[k,0]<=0): count=k k=10 else: k+=1 k=1 while (k<len(f)): B[k]= func(f,k)+2/3*h(f,k)*C[k]+1/3*h(f,k)*C[k-1] D[k]= (C[k]-C[k-1])/(3*h(f,k)) k+=1 g=A[count]+B[count]*(x-f[count,0])+C[count]*(x-f[count,0])**2+D[count]*(x-f[count,0])**3 return g def graf(f): xmin = 4.0 xmax = 4.5 dx = 0.01 xlist = mlab.frange (xmin, xmax, dx) ylist = [ispline(f,x,sweep2(f,0,0)) for x in xlist] pylab.plot (xlist, ylist) y2list = [ispline(f,x,sweep(f,0,0)) for x in xlist] pylab.plot (xlist, y2list) x = f[:,0] y = f[:,1] area = 100 scatter(x,y,s=area, marker='^', c='r') pylab.show() graf(f1)
Офлайн
А вдруг не прокатит?
Офлайн
Мало восклицательных знаков в заголовке сообщения поставил, не помогут скорее всего …
Офлайн
А задания где?
Офлайн
Прокатило.
Офлайн