SAS IML. Использование языка SAS IML. Основы IML (Лекции 2013), страница 3

метода!• Можно добавить параметр eps=1E-8, и решение будет ближе кожидаемому:call ode(r1, "FUN", c, t ,h) data="work.steps" eps=1e-8;• Но будет ли метод честно контролировать ошибку? Для этого давайтепостроим аналитическое решение, которое выражается через спец.функцию:y(x) = ( LambertW0 [c1*exp(c1-t)]+1)^(-1), где c1=1/δ -1В нашем случае δ =1E-4;• “The Lambert W function is named after J.

H. Lambert (1728 - 1777), whowas a colleague of Euler and Lagrange at the Berlin Academy of Sciences.Lambert is best known for his laws of illumination and his proof that π isirrational.“ *• … через спец. функцию, которой нет в SAS…* http://www.mathworks.com/company/newsletters/articles/stiff-differential-equations.html43Решение жесткого ОДУ. Спецфункции.• Ради методических целей разберёмся, как можно использоватьсторонние библиотеки численных методов (хотя можно было бы«упереться», и переписать вычисление LambertW0 на SAS).• Кстати, известные лектору библиотеки ч.м., которыми он сам иногдапользуется:•••http://www.netlib.org/ (много фортрана!! иногда мало документации)http://www.nr.com/ (“c”)http://www.gnu.org/software/gsl/ (доступна в любой системе, даже windows)7.23 Lambert W FunctionsLambert's W functions, W(x), are defined to be solutions of the equation W(x) \exp(W(x)) = x.

This function hasmultiple branches for x < 0; however, it has only two real-valued branches. We define W_0(x) to be the principalbranch, where W > -1 for x < 0, and W_{-1}(x) to be the other real branch, where W < -1 for x < 0. The Lambertfunctions are declared in the header file gsl_sf_lambert.h.— Function: double gsl_sf_lambert_W0 (double x)• План: подключим стороннюю библиотеку (GSL) и воспользуемсяготовой функцией44Подключение сторонних библиотек dll*• Сначала поставим уже собранную библиотеку GSL, например, отсюда:http://wiki.rglab.org/images/Data/GSL/setup_32.zip•Пропишем путь к файлу dll в переменной PATH:•Перезагрузим SAS, если он запущен.*не доступно в on-demand45Подключение сторонних библиотек• Каким образом SAS узнает о конкретной функции из конкретнойбиблиотеки?1) Сделаем специальный файл*, в котором записана сигнатура нужнойфункции:* table for sample function in gsl;* double gsl_sf_lambert_W0 (double x);ROUTINE gsl_sf_lambert_W0 MINARG=1 MAXARG=1 RETURNS=DOUBLE;ARG 1 NUM INPUT BYVALUE FORMAT=RB8.;* double gsl_sf_exp (double x);ROUTINE gsl_sf_exp MINARG=1 MAXARG=1 RETURNS=DOUBLE;ARG 1 NUM INPUT BYVALUE FORMAT=RB8.;2) Вызывающая программа:filename sascbtbl "c:\workshop\gsl_tab.sas";data work.lambert; *let's get analytical solution!;length x y arg 8;c=1E-4;set work.steps;invcm1=(1/c-1);arg=(invcm1*exp(invcm1-time));y = modulen( '*',"libgsl-0,gsl_sf_lambert_W0", arg );analyt=1/(y+1);run;* см.

help->modulen function46Подключение сторонних библиотек• При первом вызове функция modulen загружает в память указаннуюбиблиотеку.• На выходе получаем численное возвращаемое значение.• После выхода из шага data библиотека выгружается из памяти.• Можно вызывать в IML.• CAUTION:Use the correct arguments and attributes.If you use incorrect arguments or attributes, you can cause the SAS System,and possibly your operating system, to fail.• CAUTION:Using the CALL MODULE routine without a defined attribute table can causethe SAS System to fail or force you to reset your computer.You need to use an attribute table for all external functions that you want toinvoke.• Для более сложных функций (например, требующих на входструктуры/объекты, или предварительный malloc) предлагаетсяписать свою обвязку (свои dll).47Сравнение с точным решениемdata err1;set lambert;опять построим десятичныйerr=log(abs(analyt-y1))/log(10);логарифм ошибкиlabel err='log(abs(error))';run;/*plot*/goptions reset;goptions reset vsize=250pt hsize=500pt;SYMBOL1 COLOR=green INTERPOL=NONE value=x;proc gplot data=err1;y(x) = ( LambertW0 [c1*exp(c1-t)]+1)^(-1),plot err*time / vref=-8;δ =1E-4;run;quit;Здесь у экспоненты получился слишкомбольшой аргумент, ~exp(999), и мыточное решение не посчитали.где c1=1/δ -1Здесь ошибка вышла из-подконтроля численного метода.48***•На момент чтения лекции, лектору не удалось добиться от тех.

поддержки SASвнятного объяснения этого поведения функции ODE• Проблема осложняется отсутствием отладочной информации от функцииODE, или какого-либо кода возврата от неё же, который бы обращалвнимание пользователя на неконтролируемый рост ошибки.• Вроде бы можно вручную уменьшить максимальный шаг, но… он неуменьшается (вектор “h”)• Варианты решения ситуации (не полный список?):1. «Дожимать» тех. поддержку;2. Самому писать метод для интегрирования жестких ОДУ (хорошая идея!);3. Write in C:49Уравнение теплопроводности(метод конечных разностей, неявная схема)• «Экономисты любят решать уравнение теплопроводности» … (????)• Лектор не обнаружил в SAS каких-либо готовых методов для решенияуравнений в частных производных и краевых задач.• Поэтому будем писать свой! *• Формулировка задачи:* Сохин А.С., Скорик В.А.

Метод сеток решения уравнений параболического типа, Харьков, 200350Уравнение теплопроводности• Идея: введём дискретизацию по t и x, аппроксимируем производныеконечными разностями, и запишем полученную систему уравнений.• Выберем неявную схему для аппроксимации производных:Граничные условияt(j)U(i-1,j+1) U(i,j+1) U(i+1,j+1)hτU(i,j)Начальные условияПусть у нас отсутствует источник (F=0)x (i)• Сформулируем задачу в матричной форме Ax=b, где x – векторстолбец с искомыми решениями U(I,j), A – матрица, b – некоторыйвектор-столбец.51Уравнение теплопроводностидискретизованное, в матричном видеXb..U(1,1)U(1,0)+αU(0,1)A(1+2α)-α-α(1+2α)-α..U(2,1)U(2,0)-α(1+2α)..U(3,1)U(3,0)+αU(4,1)…U(1,2)αU(0,2)-1-1(1+2α)-α-α(1+2α)-α…U(2,2)0-α(1+2α)…U(3,2)αU(4,2)………-1…• Матрица А содержит только параметры разностной схемы. Столбец b содержитначальные U( *,0) и граничные U(0,*) U(4,*) условия.• В матрице много нулей.

Может не нужно их все хранить? (каждая ячейка, дажепустая – 8 байт)• Для решения такой задачи с разреженной матрицей А имеются специальныеметоды.52Разреженные матрицы•Процедура ITSOLVER получает решение разреженной линейной системы спомощью разных итеративных методов.call itsolver(x,error,it,"BICG",A,b,"MILU");• x решение системы Ax=b.• error относительная ошибка полученного решения.• iter количество итераций.Входные аргументы:method - тип итеративного метода:• "CG" - conjugate gradient algorithm.

The matrix A must be symmetric and positivedefinite.• "CGS" specifies a conjugate gradient squared algorithm, for general A.…• A – матрица в специальном формате.• b – вектор-столбец в уравнении Ax=b.• Метод для предобусловления матрицы А для более быстрой сходимости(совместимость с итеративным методом обсуждается в справке)53Разреженные матрицы/*Heat conduction equation using IML solution with finite differences */%let len=20;proc iml;h=0.1;шаги дискретизацииtau=0.1;alpha=tau/h##2;L=&len;Размер области, где получается решение = (NL+2)*hT=2;NL=int(L/h);и максимальное время = (NT+1)*tauNT=int(T/tau);Формируем разреженную матрицу A в/* valuerowcol */формате {«значение» «строка» «столбец»}do j=0 to NT-1; /*index for submatrices*/do i=1 to NL; /*fill main diagonal*/A = A //((1+2*alpha) || i+j*NL || i+j*NL);end;do i=1 to NL-1; /*fill upper&lower diagonal*/A = A //((-alpha)|| i+j*NL+1 || i+j*NL) //((-alpha)|| i+j*NL|| i+j*NL+1);end;end;do j=1 to NT-1;do i=1 to NL; /*fill 2-nd lower diagonal*/A = A //((-1)|| i+j*NL || i+(j-1)*NL);end;end;NLNL54Разреженные матрицы• Выберем начальное условие – гауссианf(x,t=0) = 1/sqrt(2*pi) *exp(-(x-L/2)^2/2)• и граничное (U(x=0)=0; U(x=L)=0).• Создадим и заполним вектор-столбец b(это обычная матрица):/* right-hand sides */b = J( NL*NT, 1, 0) ;/* add boundary conditions to b*/do i=0 to NT-1;b[i*NL+1 ,1]=b[i*NL+1,1]+alpha*0;b[i*NL+NL,1]=b[i*NL+NL,1]+alpha*0;end;/* add initial conditions to b*/do i=1 to NL-1;b[i ,1]=b[i ,1]+1/sqrt(2*3.1415926)*exp(-(L/2-i*h)##2/(2));end;• Теперь магия:call itsolver(x,error,it,"BICG",A,b,"MILU");эти настройки дают лучшую сходимость методадля несимметричной матрицы А55Уравнение теплопроводности/*break x into matrix to plot later*/SOL=J(NT*NL,3,0);do j=1 to NT; /*row*/do i=1 to NL;SOL[i+NL*(j-1),3]=x[i+NL*(j-1),1];SOL[i+NL*(j-1),1]=j*tau;SOL[i+NL*(j-1),2]=i*h;end;end;create work.temp1 from SOL;append from SOL;•Аналитическое решение для выбранного начального условия (в видегауссиана) – тоже гауссиан (известно как одномерная «предписаннаядиффузия»):•Если станд.

отклонение у гауссиана в начальном условии было σ, то:Ф(x,t)=1/ 4 + 2 2 ∗ − −02 /(4+22 )56Уравнение теплопроводностиdata analyt_temp;set temp2;retain time 1;if (col1=time) then do;analyt=1/sqrt(3.1415926*(4*time+2*1**2))*exp(-(col2-&len/2)**2/(4*time+2*1**2));output;end;run;symbol1 COLOR=green INTERPOL=join value=none;symbol2 COLOR=red INTERPOL=join value=none ;proc gplot data=analyt_temp;plot col3*col2 analyt*col2 / overlay;run;quit;57(Дискретное) Фурье-преобразование• Оно есть!• Для примера посмотрим на фильтрацию суммы синусов:%let fd=10; *discretization frequency;proc iml;t=do(-25.5,25.6,1/&fd);pi=3.1415926;y=sin(2*pi*0.2#t)+0.5*sin(2*pi*0.5#t);y1=y;do i=1 to ncol(t);y1[i]=y1[i]+RAND('NORMAL',0,0.5);end;PLT1=T(t)||T(y)||T(y1);create plot1 from plt1;append from plt1;quit;goptions reset;goptions hsize=500pt vsize=400pt;SYMBOL1 COLOR=red INTERPOL=JOIN value=none;SYMBOL2 COLOR=blue INTERPOL=JOIN value=none;proc gplot data=plot1;plot col2*col1 col3*col1/ overlay;run;quit;58Фурье-преобразованиеproc iml;use plot1;read all into plt1;NFFT=nrow(plt1); *points in signal;fft_y=fft(plt1[,2])/(NFFT /2 );*set norm for fourier transform;harm=sqrt(fft_y[,1]##2+fft_y[,2]##2);*amplitude of harmonic;fd=&fd ; *discretization frequency;f=do(0,1,2./NFFT);*frequencies for fourier-transformedsignal;Теперь можно построить:• Гармоники исходного сигнала• и зашумленного59Фурье-фильтр•Пара строчек:fil_fft_y=fft_y;fil_fft_y[loc(harm<0.1),]=0;•И обратное Фурье-преобразование:/*and inverse fourier transform*/fil_y=ifft(fil_fft_y)/2;•Теперь сбрасываем всё что нужно в наборы данных, и строим графики:fourier=T(fd/2#f)||harm||harm_1;create plot2 from fourier;append from fourier;fil_y_tab= plt1[,1]||fil_y;create plot3 from fil_y_tab;append from fil_y_tab;quit;60Фурье-фильтр• Объединим наборы данных для построения исходного, зашумленного иотфильтрованного сигнала:data merged;set plot3(rename=(col1=t1 col2=yfilt));set plot1(rename=(col1=t col2=y col3=noisy));run;• И построим график:SYMBOL1 COLOR=green INTERPOL=none value=dot;SYMBOL2 COLOR=red INTERPOL=JOIN value=none;SYMBOL3 COLOR=blue INTERPOL=JOIN value=none;proc gplot data=merged;plot noisy*t y*t yfilt*t / overlay;run;quit;61хххРазберитесь с процедурой для численного вычисления интегралов:1.

Найдите в справке, как она называется2. Напишите программу, которая вычисляет интегральный синус от +1000.Si (x )=∫0xsin(t )dtt62.