Solve PDE and Compute Partial Derivatives - MATLAB & Simulink - MathWorks Deutschland (2024)

Open Live Script

This example shows how to solve a transistor partial differential equation (PDE) and use the results to obtain partial derivatives that are part of solving a larger problem.

Consider the PDE

ut=D2ux2-DηLux.

This equation arises in transistor theory [1], and u(x,t) is a function describing the concentration of excess charge carriers (or holes) in the base of a PNP transistor. D and η are physical constants. The equation holds on the interval 0xL for times t0.

The initial condition includes a constant K and is given by

u(x,0)=KLD(1-e-η(1-x/L)η).

The problem has boundary conditions given by

u(0,t)=u(L,t)=0.

For fixed x, the solution to the equation u(x,t) describes the collapse of excess charge as t. This collapse produces a current, called the emitter discharge current, which has another constant Ip:

I(t)=[IpDKxu(x,t)]x=0.

The equation is valid for t>0 due to the inconsistency in the boundary values at x=0 for t=0 and t>0. Since the PDE has a closed-form series solution for u(x,t), you can calculate the emitter discharge current analytically as well as numerically, and compare the results.

To solve this problem in MATLAB®, you need to code the PDE equation, initial conditions, and boundary conditions, then select a suitable solution mesh before calling the solver pdepe. You either can include the required functions as local functions at the end of a file (as done here), or save them as separate, named files in a directory on the MATLAB path.

Define Physical Constants

To keep track of the physical constants, create a structure array with fields for each one. When you later define functions for the equations, the initial condition, and the boundary conditions, you can pass in this structure as an extra argument so that the functions have access to the constants.

C.L = 1;C.D = 0.1;C.eta = 10;C.K = 1;C.Ip = 1;

Code Equation

Before you can code the equation, you need to make sure that it is in the form that the pdepe solver expects:

c(x,t,u,ux)ut=x-mx(xmf(x,t,u,ux))+s(x,t,u,ux).

In this form, the PDE is

ut=x0x(x0Dux)-DηLux.

So the values of the coefficients in the equation are

  • m=0 (Cartesian coordinates with no angular symmetry)

  • c(x,t,u,ux)=1

  • f(x,t,u,ux)=Dux

  • s(x,t,u,ux)=-DηLux

Now you can create a function to code the equation. The function should have the signature [c,f,s] = transistorPDE(x,t,u,dudx,C):

  • x is the independent spatial variable.

  • t is the independent time variable.

  • u is the dependent variable being differentiated with respect to x and t.

  • dudx is the partial spatial derivative u/x.

  • C is an extra input containing the physical constants.

  • The outputs c, f, and s correspond to coefficients in the standard PDE equation form expected by pdepe.

As a result, the equation in this example can be represented by the function:

function [c,f,s] = transistorPDE(x,t,u,dudx,C)D = C.D;eta = C.eta;L = C.L;c = 1;f = D*dudx;s = -(D*eta/L)*dudx;end

(Note: All functions are included as local functions at the end of the example.)

Code Initial Condition

Next, write a function that returns the initial condition. The initial condition is applied at the first time value, and provides the value of u(x,t0) for any value of x. Use the function signature u0 = transistorIC(x,C) to write the function.

The initial condition is

u(x,0)=KLD(1-e-η(1-x/L)η).

The corresponding function is

function u0 = transistorIC(x,C)K = C.K;L = C.L;D = C.D;eta = C.eta;u0 = (K*L/D)*(1 - exp(-eta*(1 - x/L)))/eta;end

Code Boundary Conditions

Now, write a function that evaluates the boundary conditions u(0,t)=u(1,t)=0. For problems posed on the interval axb, the boundary conditions apply for all t and either x=a or x=b. The standard form for the boundary conditions expected by the solver is

p(x,t,u)+q(x,t)f(x,t,u,ux)=0.

Written in this form, the boundary conditions for this problem are

- For x=0, the equation is u+0dux=0. The coefficients are:

  • pL(x,t,u)=u,

  • qL(x,t)=0.

- Likewise for x=1, the equation is u+0dux=0. The coefficients are:

  • pR(x,t,u)=u,

  • qR(x,t)=0.

The boundary function should use the function signature [pl,ql,pr,qr] = transistorBC(xl,ul,xr,ur,t):

  • The inputs xl and ul correspond to x and u for the left boundary.

  • The inputs xr and ur correspond to x and u for the right boundary.

  • t is the independent time variable.

  • The outputs pl and ql correspond to pL(x,t,u) and qL(x,t) for the left boundary (x=0 for this problem).

  • The outputs pr and qr correspond to pR(x,t,u) and qR(x,t) for the right boundary (x=1 for this problem).

The boundary conditions in this example are represented by the function:

function [pl,ql,pr,qr] = transistorBC(xl,ul,xr,ur,t)pl = ul;ql = 0;pr = ur;qr = 0;end

Select Solution Mesh

The solution mesh defines the values of x and t where the solution is evaluated by the solver. Since the solution to this problem changes rapidly, use a relatively fine mesh of 50 spatial points in the interval 0xL and 50 time points in the interval 0t1.

x = linspace(0,C.L,50);t = linspace(0,1,50);

Solve Equation

Finally, solve the equation using the symmetry m, the PDE equation, the initial condition, the boundary conditions, and the meshes for x and t. Since pdepe expects the PDE function to use four inputs and the initial condition function to use one input, create function handles that pass in the structure of physical constants as an extra input.

m = 0;eqn = @(x,t,u,dudx) transistorPDE(x,t,u,dudx,C);ic = @(x) transistorIC(x,C);sol = pdepe(m,eqn,ic,@transistorBC,x,t);

pdepe returns the solution in a 3-D array sol, where sol(i,j,k) approximates the kth component of the solution uk evaluated at t(i) and x(j). For this problem u has only one component, but in general you can extract the kth solution component with the command u = sol(:,:,k).

u = sol(:,:,1);

Plot Solution

Create a surface plot of the solution u plotted at the selected mesh points for x and t.

surf(x,t,u)title('Numerical Solution (50 mesh points)')xlabel('Distance x')ylabel('Time t')zlabel('Solution u(x,t)')

Solve PDE and Compute Partial Derivatives- MATLAB & Simulink- MathWorks Deutschland (1)

Now, plot just x and u to get a side view of the contours in the surface plot.

plot(x,u)xlabel('Distance x')ylabel('Solution u(x,t)')title('Solution profiles at several times')

Solve PDE and Compute Partial Derivatives- MATLAB & Simulink- MathWorks Deutschland (2)

Compute Emitter Discharge Current

Using a series solution for u(x,t), the emitter discharge current can be expressed as the infinite series [1]:

I(t)=2π2Ip(1-e-ηη)n=1n2n2π2+η2/4e-dtL2(n2π2+η2/4).

Write a function to calculate the analytic solution for I(t) using 40 terms in the series. The only variable is time, but specify a second input to the function for the structure of constants.

function It = serex3(t,C) % Approximate I(t) by series expansion.Ip = C.Ip;eta = C.eta;D = C.D;L = C.L;It = 0;for n = 1:40 % Use 40 terms m = (n*pi)^2 + 0.25*eta^2; It = It + ((n*pi)^2 / m)* exp(-(D/L^2)*m*t);endIt = 2*Ip*((1 - exp(-eta))/eta)*It;end

Using the numeric solution for u(x,t) as computed by pdepe, you can also calculate the numeric approximation for I(t) at x=0 with

I(t)=[IpDKxu(x,t)]x=0.

Calculate the analytic and numeric solutions for I(t) and plot the results. Use pdeval to compute the value of u/x at x=0.

nt = length(t);I = zeros(1,nt);seriesI = zeros(1,nt);iok = 2:nt;for j = iok % At time t(j), compute du/dx at x = 0. [~,I(j)] = pdeval(m,x,u(j,:),0); seriesI(j) = serex3(t(j),C);end% Numeric solution has form I(t) = (I_p*D/K)*du(0,t)/dxI = (C.Ip*C.D/C.K)*I;plot(t(iok),I(iok),'o',t(iok),seriesI(iok))legend('From PDEPE + PDEVAL','From series')title('Emitter discharge current I(t)')xlabel('Time t')

Solve PDE and Compute Partial Derivatives- MATLAB & Simulink- MathWorks Deutschland (3)

The results match reasonably well. You can further improve the numeric result from pdepe by using a finer solution mesh.

Local Functions

Listed here are the local helper functions that the PDE solver pdepe calls to calculate the solution. Alternatively, you can save these functions as their own files in a directory on the MATLAB path.

function [c,f,s] = transistorPDE(x,t,u,dudx,C) % Equation to solveD = C.D;eta = C.eta;L = C.L;c = 1;f = D*dudx;s = -(D*eta/L)*dudx;end% ----------------------------------------------------function u0 = transistorIC(x,C) % Initial conditionK = C.K;L = C.L;D = C.D;eta = C.eta;u0 = (K*L/D)*(1 - exp(-eta*(1 - x/L)))/eta;end% ----------------------------------------------------function [pl,ql,pr,qr] = transistorBC(xl,ul,xr,ur,t) % Boundary conditionspl = ul;ql = 0;pr = ur;qr = 0;end% ----------------------------------------------------function It = serex3(t,C) % Approximate I(t) by series expansion.Ip = C.Ip;eta = C.eta;D = C.D;L = C.L;It = 0;for n = 1:40 % Use 40 terms m = (n*pi)^2 + 0.25*eta^2; It = It + ((n*pi)^2 / m)* exp(-(D/L^2)*m*t);endIt = 2*Ip*((1 - exp(-eta))/eta)*It;end% ----------------------------------------------------

References

[1] Zachmanoglou, E.C. and D.L. Thoe. Introduction to Partial Differential Equations with Applications. Dover, New York, 1986.

See Also

pdepe

Related Topics

  • Solving Partial Differential Equations
  • Solve System of PDEs

MATLAB-Befehl

Sie haben auf einen Link geklickt, der diesem MATLAB-Befehl entspricht:

 

Führen Sie den Befehl durch Eingabe in das MATLAB-Befehlsfenster aus. Webbrowser unterstützen keine MATLAB-Befehle.

Solve PDE and Compute Partial Derivatives- MATLAB & Simulink- MathWorks Deutschland (4)

Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list:

Americas

  • América Latina (Español)
  • Canada (English)
  • United States (English)

Europe

  • Belgium (English)
  • Denmark (English)
  • Deutschland (Deutsch)
  • España (Español)
  • Finland (English)
  • France (Français)
  • Ireland (English)
  • Italia (Italiano)
  • Luxembourg (English)
  • Netherlands (English)
  • Norway (English)
  • Österreich (Deutsch)
  • Portugal (English)
  • Sweden (English)
  • Switzerland
    • Deutsch
    • English
    • Français
  • United Kingdom (English)

Asia Pacific

  • Australia (English)
  • India (English)
  • New Zealand (English)
  • 中国
  • 日本 (日本語)
  • 한국 (한국어)

Contact your local office

Solve PDE and Compute Partial Derivatives
- MATLAB & Simulink
- MathWorks Deutschland (2024)
Top Articles
Latest Posts
Recommended Articles
Article information

Author: Jonah Leffler

Last Updated:

Views: 5875

Rating: 4.4 / 5 (45 voted)

Reviews: 84% of readers found this page helpful

Author information

Name: Jonah Leffler

Birthday: 1997-10-27

Address: 8987 Kieth Ports, Luettgenland, CT 54657-9808

Phone: +2611128251586

Job: Mining Supervisor

Hobby: Worldbuilding, Electronics, Amateur radio, Skiing, Cycling, Jogging, Taxidermy

Introduction: My name is Jonah Leffler, I am a determined, faithful, outstanding, inexpensive, cheerful, determined, smiling person who loves writing and wants to share my knowledge and understanding with you.