Если у вас возникли сложности с курсовой, контрольной, дипломной, рефератом, отчетом по практике, научно-исследовательской и любой другой работой - мы готовы помочь.
Предоплата всего
Подписываем
Федеральное агентство связи
ГОУ ВПО «Сибирский государственный университет
телекоммуникаций и информатики»
Уральский технический институт связи и информатики(филиал)
Курсовая работа по
Вычислительной математике
Тема: «Аппроксимация функций»
Выполнил:
Студент гр. ПЕ-71
Рябов С.А.
Екатеринбург 2008
Задание для курсовой работы.
Целью курсовой работы является комплексное применение основных вычислительных методов, изученных и апробированных на лабораторных занятиях. На первом этапе выполнения задания решается нелинейное уравнение одним из методов (по вариантам): метод половинного деления (бисекции); метод касательных; метод Вегстейна. Корень уравнения определяет интервал интегрирования дифференциального уравнения, решаемого на втором этапе задания. Здесь задача Коши решается методами интегрирования второго порядка: по средней производной и в средней точке. Полученные дискретные значения решения обрабатываются методами аппроксимации экспериментальных данных. На третьем этапе строится полином Лагранжа, интерполирующий дискретные значения на разреженной сетке узлов. На четвертом этапе сеточная функция аппроксимируется алгебраическим полиномом по методу наименьших квадратов. Соответствующая система линейных уравнений, определяющая коэффициенты полинома, решается методом квадратного корня (метод Холецкого). На заключительном этапе строятся графики сеточной функции и аппроксимирующих ее полиномов.
Решение всей задачи выполняется компьютерными методами с помощью программирования алгоритмов на одном из языков программирования (Паскаль, Бейсик). Для оценки погрешности интегрирования на втором этапе вся задача решается средствами математической системы Maple.
Постановка задачи.
Вариант № 8
(номер варианта соответствует номеру студента в списке группы)
2x^4-x^3-8
Положение корня локализовать методом прямого поиска с шагом 1. Для уточнения величины T применить метод
Вегстейна
с начальными условиями
< Y(0)=0.>2
где n – последняя цифра студенческого билета, Fs(t,) – периодическая функция с периодом =1, имеющая вид на основном отрезке периодичности [0, 1]:
Для интегрирования применить метод Эйлера второго порядка с коррекцией
< >.
Шаг интегрирования задать как где N+1 – число дискретных точек,
Построить графики сеточной функции и аппроксимирующих ее полиномов.
Решить задачу встроенными процедурами системы Maple и сравнить результаты интегрирования дифференциального уравнения.
Решение на Maple
> restart;
> with(linalg):with(plots):
pp:=(x,y)->[x,y];
Warning, the protected names norm and trace have been redefined and unprotected
Warning, the name changecoords has been redefined
> f11:= proc(t) local z ;z:=piecewise(t<0,0,t<1/2,2*t,t<=1,2-2*t,0); evalf(z);end;
> plot(f11(x),x=-1..2);
> tau:=1;
> Fs:=proc(t,tau) local x,z,y;
x:=frac(t/tau);
z:=f11(x);evalf(z);
end;
> plot(Fs(t,tau),t=-0.1..2*tau+0.1);
Полином (здесь следует задать свой вариант полинома)
> p:=2*x^4-x^3-8;
> T0:=fsolve(p,x,0..3);
> ur:=diff(Y(t),t);
Правая часть дифференциального уравнения
> Fty:=0.1*t^2-2*t*Y(t);
>
> F:=Fty+(1+2/10)*Fs(t,tau);
> RK:=dsolve({ur=F,Y(0)=0.2},Y(t),type=numeric,output=listprocedure);
> fU:=subs(RK,Y(t));
> Nt:=20;h:=T0/Nt;
> T:=array(0..Nt):U:=array(0..Nt):
> for j from 0 to Nt do
x:=j*h;z:=fU(x);T[j]:=x;U[j]:=z;
#print(x,z);
od:
> RisU:=zip(pp,T,U):
> RU:=plot(RisU):
> display(RU);
> #====================================
> RisU:=zip(pp,T,U):
> RU0:=plot(RisU,style=point,symbol=cross):
> display(RU0);
Построение интерполяционного полинома по пяти точкам
> n:=4:d:=trunc( Nt/n):X:=array(0..n):Y:=array(0..n):
> for j from 0 to n do
X[j]:=T[j*d];Y[j]:=U[j*d];
end:
> RisXY:=zip(pp,X,Y):
> GrafSym:=plot(RisXY,style=point,symbol=box,color=blue,symbolsize=16):
> display(RU0,GrafSym);
> x:='x':omega:=product((x-X[k]),k=0..n);:
> Omega:=diff(omega,x);:
> j:='j':Lagr:=sum(Y[j]*omega/(x-X[j])/(subs(x=X[j],Omega)),j=0..n);
> RLagr:=plot(evalf(Lagr),x=0..T0,color=blue):
Графики дискретных отсчетов сигнала (крестики) и аппроксимация сигнала
тригонометрическим полиномом (непрерывная линия)
> display(RLagr,RU0,GrafSym);
> Polynom:=proc(t) local z,k;global n,a;
z:=sum(a[k]*t^k,k=0..n);evalf(z);
end;
> ur:=array(0..Nt);a:=array(0..n):
> #=======Весовые функции и вид полинома=================================
> Npoint:=Nt;
Ступеньки весов
> w:='w':w:=array(0..Npoint):
> w[0]:=10:;
> for j from 1 to Npoint do
if T[j]<1 then
w[j]:=2;
else
w[j]:=1;
fi;
od:
>
>
>
> j:='j':for j from 0 to Nt do
t:=T[j];z:=Polynom(t);
ur[j]:=(z-U[j])*w[j];
end:
> UR:=convert(ur,set):
> BB:=leastsqrs(UR, {a[0],a[1],a[2],a[3],a[4]});
> a[2] = 3.530729617; a[1] = -.4472707943;a[0] = .2186165474; a[4] = 1.513446662;
a[3] = -4.387004123;
> GrafP:=plot(Polynom(x),x=0..T0):;
> display(RLagr,RU0,GrafSym,GrafP);
Формирование матрицы для метода Холецкого
> m:=n+1;A:=array(0..n,0..n):B:=array(0..n,sparse):S:=array(0..2*n,sparse);
> j:='j':k:='k':for j from 0 to Nt do
p:=1;for k from 0 to 2*n do
S[k]:=S[k]+p;p:=p*T[j];
end;
p:=1;for k from 0 to n do
B[k]:=B[k]+U[j]*p;p:=p*T[j];
end;
end:
> evalm(convert(B,vector));
> for j from 0 to n do
for k from 0 to n do
A[j,k]:=S[j+k];
end:end:
> A0:=linsolve(convert(A,matrix),convert(B,vector));
> for j from 0 to n do a[j]:=A0[j+1];end:
> t:='t':GrafG:=plot(Polynom(t),t=0..T0,color=brown):;
> display(RLagr,RU0,GrafSym,GrafG,GrafP);
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++=
> fn:=`C:\\work\\rezult.txt`;
> L:=readdata(fn,3):;
Nstrok:=vectdim(L);
> U_n:=array(1..Nstrok);:
T_n:=array(1..Nstrok);
> for j from 1 to Nstrok do
T_n[j]:=L[j,1];
U_n[j]:=L[j,2];
#print(j,T_n[j],U_n[j]);
od:
> u1:=zip(pp,T_n,U_n):
> RU1:=plot(u1,style=point,symbol=cross,color=black):
> display(RU,RU1);
Сравнение результатов на одинаковых сетках
>
> printf("%s",` № t U_map U_pas разн`);
for k from 0 to Nt do t:=T[k]:del:=U[k]-U_n[k+1];
printf("% 3.0f % 6.2f % 8.4f % 8.4f % 8.4f \n",k,t,U[k],U_n[k+1],del):
end:;
№ t U_map U_pas разн
0 0.00 .2000 .2000 0.0000
1 .08 .2061 .2182 -.0121
2 .16 .2241 .2479 -.0238
3 .23 .2536 .2882 -.0346
4 .31 .2935 .3377 -.0442
5 .39 .3426 .3949 -.0523
6 .47 .3993 .4556 -.0563
7 .55 .4572 .4924 -.0352
8 .62 .4941 .5068 -.0127
9 .70 .5087 .5010 .0077
10 .78 .5031 .4774 .0257
11 .86 .4796 .4387 .0409
12 .93 .4410 .3879 .0531
13 1.01 .3905 .3468 .0437
14 1.09 .3485 .3228 .0257
15 1.17 .3238 .3137 .0101
16 1.25 .3142 .3174 -.0032
17 1.32 .3173 .3314 -.0141
18 1.40 .3310 .3536 -.0226
19 1.48 .3530 .3748 -.0218
20 1.56 .3736 .3755 -.0019
Решение на Pascal
program kurs;
uses crt;
const
nt=20;
u0=0.2;
type
massiv= array[0..nt] of real;
var h,t0,z,root: real;
fu: text;
u: massiv;
function pol(t:real):real;
begin
pol:=2*t*t*t*t-t*t*t-8;
end;
{*********}
function deriv(t:real):real;
begin
deriv:=8*t*t*t-3*t*t;
end;
{*********}
procedure vegst(a,b:real; var x:real);
var
x0,x1,y1,y0,y2,x2:real;
delta:real;
i:integer;
begin
x0:=1.5;
y0:=x0;
x1:=2.2;
y1:=x1;
i:=1;
while i<>5 do
begin
y2:=x1+pol(x1);
x2:=y2-(y2-y1)*(y2-x1)/((y2-y1)-(x1-x0));
writeln('x1= ',x1:10:3,' x2= ',x2:10:3);
y1:=y2;
y2:=x2+pol(x2);
x0:=x1;
x1:=x2;
i:=i+1;
end;
x:=x1;
end;
{**************}
function signal(e: real): real;
begin
if e<0 then
z:=0
else
if e<1/2 then
z:=2*e
else
if e<=1 then
z:=2-2*e
else
z:=0;
signal:=z;
end;
function Fs(t: real): real;
begin
z:=t-trunc(t);
Fs:=signal(z);
end;
function RHS(t,y : real): real;
begin
RHS:= 0.1*t*t-2*t*y + (1+2/10)*Fs(t);
end;
{**********}
procedure difur(var u: massiv);
var
t: real;
k: integer;
begin
u[0]:=u0;
writeln(fu,0.0:8:2,u0:10:4);
for k:=1 to nt do
begin
t:=k*h;
u[k]:=u[k-1]+h*rhs(t+h/2,u[k-1]+h/2*rhs(t,u[k-1]));
writeln(fu,t:8:2,u[k]:10:4);
end;
end;
begin
clrscr;
vegst(0,3,root);
t0:=root;
h:=t0/nt;
assign(fu,'c:\work\rezul.txt');
rewrite(fu);
difur(u);
close(fu);
repeat until keypressed;
end.
Перечень литературы