Применение численных методов для решения уравнений с частными производными

/>

TABL=TABL'

По данным таблицы запишем и решим нормальную систему МНК-метода:

2) дл многочлена второй степени

S2=[n TABL(7,1) TABL(7,3);TABL(7,1) TABL(7,3) TABL(7,6);TABL(7,3) TABL(7,6) TABL(7,7)] матрица коэффициентов

T2=[TABL(7,2);TABL(7,4);TABL(7,5)] вектор правых частей coef2=S2T2 решение нормальной системы МНК

A2=coef2(3);B2=coef2(2);C2=coef2(1); коэффициенты многочлена 2-ой степени

Дл построения графиков функции y2=A2*x^2+B2*x+C2 с найденными коэффициентами зададим вспомогательный вектор абсциссы xi, а затем вычислим элементы векторов g1=A1*xi+B1 и g2=A2*xi^2+B2*xi+C2 :

h=0.05;

xi=min(X):h:max(X);

g2=log(1./(A2*xi.^2+B2*xi+C2));

plot(X,Y_o,'*k',xi,g2);grid

coef2=polyfit(X,Y,2) коэффициенты многочлена второй степени

Для построения графиков зададим вспомогательный вектор абсциссы xi, а затем c помощью функции polyval вычислим элементы векторов g1 и g2:

pause;

xi=min(X):0.1:max(X);

g2=polyval(coef2,xi);

plot(X,Y_o,'*k',xi,log(1./g2));grid

Очевидно, что построенные таким способом графики совпадут с полученными ранее.

Для того, чтобы определить величину среднеквадратичного уклонения, вычислим суммы квадратов уклонений g1(x) и g2(x) от таблично заданной функции в узлах таблицы X а затем

G2=polyval(coef2,X);

delt2=sum((Y-G2).^2); delt2=sqrt(delt2/5)

Последние две строки можно заменить другими, если использовать функцию mean , вычисляющую среднее значение:

delt2=mean(sum((Y-G2).^2))

Y = 0.0378 0.0653 0.3789 1.0353 0.5172 0.9765

Y_o = 0.0378 0.0653 0.3789 1.0353 0.5172 0.9765

Y = 0.9629 0.9368 0.6846 0.3551 0.5962 0.3766

n = 6

TABL =

0 0.9629 0 0 0 0 0

0.5000 0.9368 0.2500 0.4684 0.2342 0.1250 0.0625

1.0000 0.6846 1.0000 0.6846 0.6846 1.0000 1.0000

1.5000 0.3551 2.2500 0.5327 0.7990 3.3750 5.0625

2.0000 0.5962 4.0000 1.1924 2.3848 8.0000 16.0000

2.5000 0.3766 6.2500 0.9416 2.3539 15.6250 39.0625

7.5000 3.9122 13.7500 3.8196 6.4565 28.1250 61.1875

S2 =

6.0000 7.5000 13.7500

7.5000 13.7500 28.1250

13.7500 28.1250 61.1875

T2 =

3.9122

3.8196

6.4565

coef2 =

1.0178

-0.4243

0.0718

coef2 =

0.0718 -0.4243 1.0178

delt2 =

0.1199

delt2 =

0.0719

Численные методы решения задачи Коши для обыкновенных дифференциальных уравнений

Эйлера явная

function dy=func(x,y)

dy=2*x*y

clear

X=[0.00000 0.10000 0.20000 0.30000 0.40000 0.50000];

Y=exp((X).^2);

Y0=input('Значение функции в точке 0 = ');

Y_n1=Y0;

for n=1:length(X)-1

Y_n1=Y_n1+0.1*2*X(n)*Y_n1;

Y_n(n)=Y_n1;

end

X1=0.00000:0.01:0.50000;

Y_sot=Y0;

for n=1:length(X1)

Y_sot=Y_sot+0.01*2*X1(n)*Y_sot;

Y_sto(n)=Y_sot;

end

X2=0.00000:0.001:0.50000;

Y_tys=Y0;

for n=1:length(X2)

Y_tys=Y_tys+0.001*2*X2(n)*Y_tys;

Y_ts(n)=Y_tys;

end

disp(' X Y h=0.1 h=0.01 h=0.001')

G1=Y_sto(10:10:end);

TABL=[X;Y;Y0,Y_n;...

Y_sto(1),G1;...

Y_ts(1),Y_ts(100:100:end)];

TABL=TABL';disp(TABL)

Значение функции в точке 0 = 1

X Y h=0.1 h=0.01 h=0.001

0 1.0000 1.0000 1.0000 1.0000

0.1000 1.0101 1.0000 1.0090 1.0099

0.2000 1.0408 1.0200 1.0387 1.0406

0.3000 1.0942 1.0608 1.0907 1.0938

0.4000 1.1735 1.1244 1.1683 1.1730

0.5000 1.2840 1.2144 1.2766 1.2833

Симметричная

clear

X=[0.00000 0.10000 0.20000 0.30000 0.40000 0.50000];

Y=exp((X).^2);

Y0=input('Значение функции в точке 0 = ');

Y_n1=Y0;

for n=1:length(X)-1

Y_n1=Y_n1*(1+0.1*X(n))/(1-0.1*(X(n)+0.1));

Y_n(n)=Y_n1;

end

X1=0.00000:0.01:0.50000;

Y_sot=Y0;

for n=1:length(X1)-1

Y_sot=Y_sot*(1+0.01*X1(n))/(1-0.01*(X1(n)+0.01));

Y_sto(n)=Y_sot;

end

X2=0.00000:0.001:0.50000;

Y_tys=Y0;

for n=1:length(X2)-1

Y_tys=Y_tys*(1+0.001*X2(n))/(1-0.001*(X2(n)+0.001));

Y_ts(n)=Y_tys;

end

disp(' X Y h=0.1 h=0.01 h=0.001')

G1=Y_sto(10:10:end);

TABL=[X;Y;Y0,Y_n;...

Y_sto(1),G1;...

Y_ts(1),Y_ts(100:100:end)];

TABL=TABL';disp(TABL)

Значение функции в точке 0 = 1

X Y h=0.1 h=0.01 h=0.001

0 1.0000 1.0000 1.0001 1.0000

0.1000 1.0101 1.0101 1.0101 1.0101

0.2000 1.0408 1.0410 1.0408 1.0408

0.3000 1.0942 1.0947 1.0942 1.0942

0.4000 1.1735 1.1745 1.1735 1.1735

0.5000 1.2840 1.2858 1.2840 1.2840

Эйлера неявная

clc

clear all

h_1=0.1;

X=0:h_1:0.5;

Y=exp(X.^2);

Yn=Y(1);

Y2=Yn+h_1*2*X(1);

или Y2=input('Введите значение Yn в точке X=0 =')

y_1(1)=Y2;

for i = 1:(length(Y)-1)

y_1(i+1)=y_1(i)/(1-h_1*2*X(i+1));

end

h_2=0.01;

X_2=0:h_2:0.5;

Y_2=exp(X_2.^2);

Y2=Yn+h_2*2*X(1);

y_2(1)=Y2;

for i = 1:(length(Y_2)-1)

y_2(i+1)=y_2(i)/(1-h_2*2*X_2(i+1));

end

h_3=0.001;

X_3=0:h_3:0.5;

Y_3=exp(X_3.^2);

Y2=Yn+h_3*2*X(1);

y_3(1)=Y2;

for i = 1:(length(Y_3)-1)

y_3(i+1)=y_3(i)/(1-h_3*2*X_3(i+1));

end

for k=1:5

r1(k)=y_2(k*10+1);

r2(k)=y_3(k*100+1);

end

TABL=[X; Y;... ... означает перенос строки

y_1;...

y_2(1),r1;...

y_3(1),r2];

TABL=TABL'

plot(X,Y,'-',X,y_1,X,[y_2(1),r1],X,[y_3(1),r2])

f=ode23('y_1',[0 0.5],1)

TABL =

0 1.0000 1.0000 1.0000 1.0000

0.1000 1.0101 1.0204 1.0111 1.0102

0.2000 1.0408 1.0629 1.0430 1.0410

0.3000 1.0942 1.1308 1.0977 1.0945

0.4000 1.1735 1.2291 1.1787 1.1740

0.5000 1.2840 1.3657 1.2916 1.2848

Задача
Коши

function [ output_args ] = koshi( input_args )

KOSHI Summary of this function goes here

Detailed explanation goes here

tspan=[0,1];

y0=[1.421,1];

[t,y]=ode45(@F,tspan,y0);

ode45(@F,tspan,y0);

hold on

plot([0 1],[1 1])

Подбор Альфа методом секущих

a=1;

y0=[1,a];

tspan=[0,1];

psi_old=a-1;

a_old=0.5;

i = 1;

eps = 1;

while (eps >= 0.000001) & (i < 10000)

[t,y]=ode45(@F,tspan,y0);

ode45(@F,tspan,y0)

psi=y(end,2)-1;

a_new=a-psi*(a-a_old)/(psi-psi_old)

eps = abs(psi);

a_old = a;

a = a_new;

y0=[1,a];

psi_old = psi

i = i + 1;

end

i

a_new = 0.5000

psi_old = -0.3935

a_new = 1.4655

psi_old = -0.8161

a_new = 1.4184

psi_old = 0.0419

a_new = 1.4215

psi_old = -0.0030

a_new = 1.4215

psi_old = -4.1359e-006

a_new = 1.4215

psi_old = 4.2046e-010

i = 7

Генерация матрицы 10*10 из элементов равномерно распределённых 1..10

function [ output_args ] = ravnomern10_10_1_10( input_args )

RAVNOMERN10_10_1_10 Summary of this function goes here

Detailed explanation goes here

round(rand(10,10)*9+1)

8 2 7 7 5 3 8 9 4 2

9 10 1 1 4 7 3 3 8 1

2 10 9 3 8 7 6 8 6 6

9 5 9 1 8 2 7 3 6 8

7 8 7 2 3 2 9 9 9 9

2 2 8 8 5 5 10 4 4 2

4 5 8 7 5 10 6 3 8 6

6 9 5 4 7 4 2 3 8 5

10 8 7 10 7 6 2 7 4 1

10 10 3 1 8 3 3 5 6 4

Решение краевой задачи методом сеток для УЧП.

e=0.01;

h=sqrt(e);

x=0:h:1;

y=0:h:1;

v=ones(11,11);

v(1,:)=0;

v(end,:)=1;

v(:,1)=(0:h:1)';

v(:,end)=(0:h:1)';

v=v.*((1*9+sum(0:h:1)+sum(0:h:1))/40)

v(1,:)=0;

v(end,:)=1;

v(:,1)=(0:h:1)';

v(:,end)=(0:h:1)';

surf(v);

d = e+1;

i=1

while d > e & i < 100

v1=v;

v1(1:10,:)=v1(2:11,:);

v1(11,:)=v(1,:);

v2=v;

v2(2:11,:)=v2(1:10,:);

v2(1,:)=v(11,:);

v3=v;

v3(:,1:10)=v3(:,2:11);

v3(:,11)=v(:,1);

v4=v;

v4(:,2:11)=v4(:,1:10);

v4(:,1)=v(:,11);

v_new=(v1+v2+v3+v4)/4;

d = max(max(abs(v(2:end-1,2:end-1)-v_new(2:end-1,2:end-1))))

v=v_new;

v(1,:)=0;

v(end,:)=1;

v(:,1)=(0:h:1)';

v(:,end)=(0:h:1)';

pause(0.5);

surf(v);

i = i + 1;

end;

Результат работы программы:

i = 1

d = 0.2250

d = 0.0875

d = 0.0500

d = 0.0344

d = 0.0297

d = 0.0245

d = 0.0199

d = 0.0175

d = 0.0154

d = 0.0137

d = 0.0120

d = 0.0108

d = 0.0093

HELM ВЫЧИСЛЯЕТ МЕТОДОМ МОНТЕ-КАРЛО (АЛГОРИТМ "БЛУЖДАНИЙ ПО СЕТКЕ") РЕШЕНИЕ ЗАДАЧИ ДИРИХЛЕ ДЛЯ УРАВНЕНИЯ ГЕЛЬМГОЛЬЦА В ЗАДАННОЙ ТОЧКЕ (x,y) ПРЯМОУГОЛЬНОЙ ОБЛАСТИ D, ОПРЕДЕЛЕННОЙ ГРАНИЦАМИ.

Код программ:

ЗАПИС КРАЕВЫХ УСЛОВИЙ В ЗАДАЧЕ

ДИРИХЛЕ ДЛЯ УРАВНЕНИЯ ГЕЛЬМГОЛЬЦА

function yp=funch(x,y);

if x=0,yp=y;end;

if y=0,yp=0;end;

if y=1,yp=1;end;

if x=1,yp=y;end;

function [z1,z2,z3]=helm(c,fun,xm,ym,gr,x0,y0,h,n);

HELM ВЫЧИСЛЯЕТ МЕТОДОМ МОНТЕ-КАРЛО (АЛГОРИТМ "БЛУЖДАНИЙ ПО СЕТКЕ")

РЕШЕНИЕ ЗАДАЧИ ДИРИХЛЕ ДЛЯ УРАВНЕНИЯ ГЕЛЬМГОЛЬЦА В ЗАДАННОЙ

ТОЧКЕ (x,y) ПРЯМОУГОЛЬНОЙ ОБЛАСТИ D, ОПРЕДЕЛЕННОЙ ГРАНИЦАМИ

0<=x<=xm, 0<=y<=ym

(УЧП) Uxx+Uyy-c*U=F(x,y)

(ГУ) U|г=G(x,y)

Входные данные:

c - КОЭФФИЦИЕНТ (функциональный) ЛЕВОЙ ЧАСТИ УЧП;

fun - ФУНКЦИЯ F(x,y) В ПРАВОЙ ЧАСТИ УЧП (ФУНКЦИЯ ПОЛЬЗОВАТЕЛЯ);

xm,ym - ГРАНИЦЫ ПРЯМОУГОЛЬНОЙ ОБЛАСТИ;

gr - ГРАНИЧНЫЕ УСЛОВИЯ (ФУНКЦИЯ ПОЛЬЗОВАТЕЛЯ);

x0,y0 - КООРДИНАТЫ ТОЧКИ, В КОТОРОЙ ИЩЕТСЯ РЕШЕНИЕ;

h - ШАГ СЕТКИ (ЗАДАЕТСЯ ПОЛЬЗОВАТЕЛЕМ);

n - ЧИСЛО ТРАЕКТОРИЙ.

Выходные данные:

z1 - ПРИБЛИЖЕННОЕ ЗНАЧЕНИЕ РЕШЕНИЯ;

z2 - ВЕРОЯТНАЯ ОШИБКА;

z3 - ВЕРХНЯЯ ГРАНИЦА ОШИБКИ.

Обращение: [z1,z2,z3]=helm(c,fun,xm,ym,gr,x0,y0,h,n)

global z

z=[];

i0=round(x0/h);

j0=round(y0/h);

im=round(xm/h);

jm=round(ym/h);

disp(' ')

disp(' Подождите, идет расчет.')

for count=1:n,

x=x0;y=y0;

i=i0;j=j0;

q=1;

tmp=4+eval(c)*h^2;

s=h^2*eval(fun)/tmp;

while all([i,j,im-i,jm-j]),

p=[0,1/4];p=[p,p(2)];

p=[p,1/4]; p=[p,p(4)];

alf=rand;

pp=max(find(alf>cumsum(p)));

if pp==1,j=j+1;end

if pp==2,j=j-1;end

if pp==3,i=i+1;end

if pp==4,i=i-1;end

x=i*h;y=j*h;

q=q*4/tmp;

s=s+q*h^2*eval(fun)/tmp;

end

s=s+q*feval(gr,x,y);

z=[z,s];

end

disp(' ');

disp(' РЕШЕНИЕ ЗАДАЧИ:');

disp(' ============================= ');

disp(' ')

disp(' при числе траекторий');disp(n);

disp('значение в точке с координатами ');

disp(' x0 y0');

disp([x0,y0]);

z1=mean(z);disp(' ПРИБЛИЖЕННОГО РЕШЕНИЯ - ');disp(z1);

z2=0.6745*std(z)/sqrt(n);disp(' ВЕРОЯТНОЙ ОШИБКИ - ');disp(z2);

z3=z2*3/0.6745;disp(' ВЕРХНЕЙ ГРАНИЦЫ ОШИБКИ - ');disp(z3);

ОБРАЩЕНИЯ К ФУНКЦИИ HELM:

global z

c='1';

f='0';

xm=1;ym=1;

gr='funch';

x0=0.6;y0=0.7;

h=0.1;

n=600;

[z1,z2,z3]=helm(c,f,xm,ym,gr,x0,y0,h,n);

Результат работы программы:

РЕШЕНИЕ ЗАДАЧИ:

при числе траекторий 600

значение в точке с координатами x0 y0 0.6000 0.7000

ПРИБЛИЖЕННОГО РЕШЕНИЯ - 0.2958

ВЕРОЯТНОЙ ОШИБКИ - 0.0089

ВЕРХНЕЙ ГРАНИЦЫ ОШИБКИ - 0.0397