Создание символьного PDE, представляющего лапласиан переменной u(x,y).

syms u(x,y)
pdeeq = laplacian(u,[x y]) == -3

pdeeq(x, y) = 
$$\frac{{\partial }^2}{\partial x^2}\mathrm{ }u\left(x,y\right)+\frac{{\partial }^2}{\partial y^2}\mathrm{ }u\left(x,y\right)=-3$$

Извлеките коэффициенты PDE.

coeffs = pdeCoefficients(pdeeq,u)

coeffs = struct with fields:
    a: 0
    c: [4x1 double]
    m: 0
    d: 0
    f: -3

pdeCoefficients преобразует символьный PDE в скалярное уравнение PDE вида

и извлекает коэффициенты m, d, c, a, и f в структуру coeffs. Для 2-D систем $$N$$ уравнений, c - тензор с $${}^{\mathrm{4N2}}$$ элементами. Для получения дополнительной информации см. c Коэффициент для specifyCoefficients (Частичный Отличительный Комплект инструментов Уравнения). В этом случае $$N\mathrm{=}$$1, поэтому коэффициентc является вектором строки 4 на 1, который представляет c$${}_{\mathrm{xx}}$$, c$${}_{\mathrm{xy}}$$, c$${}_{\mathrm{yx}}$$, и c$${}_{\mathrm{yy}}$$.

coeffs.c

ans = 4×1

    -1
     0
     0
    -1

Решите 2-D однородное уравнение Лапласа в области, ограниченной единичной окружностью. Используйте pdeCoefficients функция для извлечения коэффициентов уравнения Лапласа.

Создайте модель PDE и включите геометрию.

model = createpde;
geometryFromEdges(model,@circleg);

Постройте график геометрии и отобразите метки кромок для использования в определении граничного условия.

figure
pdegplot(model,'EdgeLabels','on')
axis equal

Создание символического выражения pdeeq это представляет уравнение Лапласа.

syms u(x,y)
pdeeq = laplacian(u,[x y])

pdeeq(x, y) = 
$$\frac{{\partial }^2}{\partial x^2}\mathrm{ }u\left(x,y\right)+\frac{{\partial }^2}{\partial y^2}\mathrm{ }u\left(x,y\right)$$

Извлеките коэффициенты уравнения Лапласа.

coeffs = pdeCoefficients(pdeeq,u)

coeffs = struct with fields:
    a: 0
    c: [4x1 double]
    m: 0
    d: 0
    f: 0

Укажите коэффициенты модели PDE.

specifyCoefficients(model,'m',coeffs.m,'d',coeffs.d, ...
    'c',coeffs.c,'a',coeffs.a,'f',coeffs.f);

Примените граничное условие Дирихле $$u(x,y{)}^{\mathrm{}}={}^{\mathrm{}}$$x4 + y4 ко всем ребрам 4 которые образуют окружность.

applyBoundaryCondition(model,'dirichlet','Edge',1:4,'u',@(region,state) region.x.^4 + region.y.^4);

Создайте сетку по умолчанию для геометрии.

generateMesh(model);

Решите PDE и постройте график решения.

results = solvepde(model);
pdeplot(model,'XYData',results.NodalSolution)

Решите уравнение 2-го Пуассона в области, ограниченной кругом единицы. Форма расхождения PDE имеет непостоянную f коэффициент. Можно решить PDE, извлекая коэффициенты с помощью pdeCoefficients и определение коэффициентов в модели PDE с использованием specifyCoefficients.

Создайте модель PDE и включите геометрию.

model = createpde;
geometryFromEdges(model,@circleg);

Создание символического выражения pdeeq это уравнение Пуассона.

syms u(x,y)
pdeeq = diff(u,x,x) + diff(u,y,y) - 1/(x^2 + y^2)

pdeeq(x, y) = 
$$\frac{{\partial }^2}{\partial x^2}\mathrm{ }u\left(x,y\right)-\frac{1}{x^2+y^2}+\frac{{\partial }^2}{\partial y^2}\mathrm{ }u\left(x,y\right)$$

Извлеките коэффициенты уравнения.

coeffs = pdeCoefficients(pdeeq,u)

coeffs = struct with fields:
    a: 0
    c: [4x1 double]
    m: 0
    d: 0
    f: @makeCoefficient/coefficientFunction

f коэффициент не является константой. Он отображается как @makeCoefficient/coefficientFunction, что указывает на то, что он был возвращен из рабочей области какой-либо внутренней функции. Чтобы показать формулу за дескриптором функции, используйте coeffs.f('show').

coeffs.f('show')

ans = 
$$\frac{1}{x^2+y^2}$$

Укажите коэффициенты модели PDE.

specifyCoefficients(model,'m',coeffs.m,'d',coeffs.d, ...
    'c',coeffs.c,'a',coeffs.a,'f',coeffs.f);

Примените граничное условие Дирихле $$u(x,y\mathrm{)}$$ = 0 на всех 4 кромках, образующих окружность.

applyBoundaryCondition(model,'dirichlet','Edge',1:4,'u',0);

Создайте сетку по умолчанию для геометрии.

generateMesh(model);

Решите проблему PDE. Постройте график узлового решения с помощью опции 'XYData' в pdeplot функция.

results = solvepde(model);
pdeplot(model,'XYData',results.NodalSolution)

Постройте график градиента решения в узловых местах с помощью опции 'FlowData'.

pdeplot(model,'FlowData',[results.XGradients results.YGradients])

Создайте PDE без формы расхождения.

syms u(x,y)
pdeeq = diff(u,x,x) + cos(x+y)/4*diff(u,x,y) + 1/2*diff(u,y,y)

pdeeq(x, y) = 
$$\frac{{\partial }^2}{\partial x^2}\mathrm{ }u\left(x,y\right)+\frac{\cos \left(x+y\right)\hspace{0.17em}\frac{\partial }{\partial y}\mathrm{ }\frac{\partial }{\partial x}\mathrm{ }u\left(x,y\right)}{4}+\frac{\frac{{\partial }^2}{\partial y^2}\mathrm{ }u\left(x,y\right)}{2}$$

Извлеките его коэффициенты. Когда pdeCoefficients не удается преобразовать PDE в форму расхождения

$$\frac{{}^{\mathrm{}}}{{}^{\mathrm{}}}\frac{}{}m\partial 2u\partial t2+d\partial u\partial t-\nabla \cdot (c\otimes \nabla u)+$$au = f,

он выдает предупреждающее сообщение, но по-прежнему выдает выходные данные.

coeffs = pdeCoefficients(pdeeq,u)

Warning: After extracting m, d, and c, some gradients remain. Writing all remaining terms to f.

coeffs = struct with fields:
    a: 0
    c: @makeCoefficient/coefficientFunction
    m: 0
    d: 0
    f: @makeCoefficient/coefficientFunction

Отображение дескрипторов функции в извлеченных коэффициентах c и f, используйте опцию 'show'. Все оставшиеся градиенты в PDE записываются в f коэффициент.

coeffs.c('show')

ans = 
$$\left(\begin{array}{c}-1\\ \frac{1}{8}-\frac{\cos \left(x+y\right)}{4}\\ \frac{1}{8}-\frac{\cos \left(x+y\right)}{4}\\ -\frac{1}{2}\end{array}\right)$$

coeffs.f('show')

ans = 
$$\frac{\sin \left(x+y\right)\hspace{0.17em}\frac{\partial }{\partial x}\mathrm{ }u\left(x,y\right)}{4}+\frac{\sin \left(x+y\right)\hspace{0.17em}\frac{\partial }{\partial y}\mathrm{ }u\left(x,y\right)}{4}-\frac{\cos \left(x+y\right)\hspace{0.17em}\frac{\partial }{\partial y}\mathrm{ }\frac{\partial }{\partial x}\mathrm{ }u\left(x,y\right)}{4}+\frac{\frac{\partial }{\partial y}\mathrm{ }\frac{\partial }{\partial x}\mathrm{ }u\left(x,y\right)}{4}$$

Поскольку PDE не имеет формы расхождения, требуемой для панели инструментов PDE, панель инструментов не сможет решить эту PDE.

Решите зависящее от времени волновое уравнение в области, ограниченной единичной окружностью. Волновое уравнение представляет собой PDE со скалярной функцией $$u(t,x,$$ y), которая зависит от $$$$времени t и $$$$координат x $$$$и y.

Создайте модель PDE и включите геометрию.

model = createpde(1);
geometryFromEdges(model,@circleg);

Создайте символьный PDE, представляющий волновое уравнение.

syms u(t,x,y)
pdeeq = diff(u,t,t) - laplacian(u,[x y])

pdeeq(t, x, y) = 
$$\frac{{\partial }^2}{\partial t^2}\mathrm{ }u\left(t,x,y\right)-\frac{{\partial }^2}{\partial x^2}\mathrm{ }u\left(t,x,y\right)-\frac{{\partial }^2}{\partial y^2}\mathrm{ }u\left(t,x,y\right)$$

Извлеките коэффициенты PDE.

coeffs = pdeCoefficients(pdeeq,u)

coeffs = struct with fields:
    a: 0
    c: [4x1 double]
    m: 1
    d: 0
    f: 0

Укажите коэффициенты модели PDE.

specifyCoefficients(model,'m',coeffs.m,'d',coeffs.d, ...
    'c',coeffs.c,'a',coeffs.a,'f',coeffs.f);

Задайте начальные условия зависящей от времени задачи для всей геометрии.

setInitialConditions(model,0,1);

Примените граничное условие Дирихле $$u(x,y\mathrm{)}$$ = 0 на всех 4 кромках, образующих окружность.

applyBoundaryCondition(model,'dirichlet','Edge',1:4,'u',0);

Создайте сетку по умолчанию для геометрии.

generateMesh(model);

Найдите решения зависящего от времени PDE в диапазоне времени от 0 с до 50 с интервалом 2 с.

results = solvepde(model,linspace(0,2,50));

Постройте график решения волнового уравнения для каждого интервала 5 с.

figure;
for k = 1:10
    subplot(5,2,k);
    pdeplot(model,'XYData',results.NodalSolution(:,k*5))
    axis equal
end  

Решите систему из двух PDE второго порядка. Вы можете решить систему PDE, извлекая коэффициенты PDE символически с помощью pdeCoefficients, преобразование коэффициентов в числа с двойной точностью с помощью pdeCoefficientsToDoubleи определение коэффициентов в модели PDE с использованием specifyCoefficients.

Система ПДЭ представляет собой прогиб зажатой конструкционной пластины при равномерной нагрузке давлением. Система PDE с зависимыми переменными $${}_{\mathrm{u1}}$$ и $${}_{\mathrm{u2}}$$ задается

$$-{}^{\mathrm{}}{}_{\mathrm{}}{}_{\mathrm{}}\mathrm{}$$∇2u1+u2=0,

$$-{}^{\mathrm{}}{}_{\mathrm{}}$$D∇2u2=p,

где $$D$$ - жесткость пластины при изгибе,

$$D\frac{={}^{\mathrm{}}}{\mathrm{}\mathrm{}\mathrm{Eh312}{(}^{\mathrm{}}}$$1-start2),

и $$E$$ - модуль упругости, $$\mathrm{start}$$- отношение Пуассона, $$h$$ - толщина пластины, $${}_{\mathrm{u1}}$$ - поперечное отклонение пластины, $$p$$ - нагрузка давления.

Создайте модель PDE для системы двух уравнений.

model = createpde(2);

Создайте квадратную геометрию. Укажите длину стороны квадрата. Затем включите геометрию в модель PDE.

len = 10.0;
gdm = [3 4 0 len len 0 0 0 len len]';
g = decsg(gdm,'S1',('S1')');
geometryFromEdges(model,g);

Укажите значения физических параметров системы. Пусть внешнее давление $$p$$ является символьной переменной pres это может принимать любое значение.

E = 1.0e6;
h_thick = 0.1;
nu = 0.3;
D = E*h_thick^3/(12*(1 - nu^2));
syms pres

Объявите систему PDE как систему символьных уравнений. Извлеките коэффициенты PDE и верните их в символической форме.

syms u1(x,y) u2(x,y)
pdeeq = [-laplacian(u1) + u2; -D*laplacian(u2) - pres];
symCoeffs = pdeCoefficients(pdeeq,[u1 u2],'Symbolic',true)

symCoeffs = struct with fields:
    m: [1x1 sym]
    a: [2x2 sym]
    c: [4x4 sym]
    f: [2x1 sym]
    d: [1x1 sym]

Отображение коэффициентов m, a, c, f, и d.

structfun(@disp,symCoeffs)

$$0$$

$$\left(\begin{array}{cc}0& 1\\ 0& 0\end{array}\right)$$

$$\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& \frac{25000}{273}& 0\\ 0& 0& 0& \frac{25000}{273}\end{array}\right)$$

$$\left(\begin{array}{c}0\\ \mathrm{pres}\end{array}\right)$$

$$0$$

Заменить значение на pres с использованием subs функция. С момента выхода subs являются символическими объектами, используйте pdeCoefficientsToDouble для преобразования коэффициентов в double тип данных, который делает их действительными входными данными для панели инструментов PDE.

symCoeffs = subs(symCoeffs,pres,2);
coeffs = pdeCoefficientsToDouble(symCoeffs)

coeffs = struct with fields:
    a: [4x1 double]
    c: [16x1 double]
    m: 0
    d: 0
    f: [2x1 double]

Укажите коэффициенты PDE для модели PDE.

specifyCoefficients(model,'m',coeffs.m,'d',coeffs.d, ...
    'c',coeffs.c,'a',coeffs.a,'f',coeffs.f);

Задайте жесткость пружины. Задание граничных условий путем определения распределенных пружин на всех четырех кромках.

k = 1e7;
bOuter = applyBoundaryCondition(model,'neumann','Edge',(1:4), ...
    'g',[0 0],'q',[0 0; k 0]);

Укажите размер сетки геометрии и создайте сетку для модели PDE.

hmax = len/20;
generateMesh(model,'Hmax',hmax);

Решите модель.

res = solvepde(model);

Доступ к решению в узловых местоположениях.

u = res.NodalSolution;

Постройте график поперечного отклонения пластины.

numNodes = size(model.Mesh.Nodes,2);
figure;
pdeplot(model,'XYData',u(1:numNodes),'contour','on')
title 'Transverse Deflection'

Найдите поперечное отклонение в центре пластины.

wMax = min(u(1:numNodes,1))

wMax = -0.2763

Сравните результат с отклонением в центре пластины, вычисленным аналитически.

pres = 2;
wMax = -.0138*pres*len^4/(E*h_thick^3)

wMax = -0.2760